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Abstract 

We  develop  divergence-conforming  B-splinc  discretizations  for  the  numerical  solution 
of  the  Darcy-Stokes-Brinkman  equations.  These  discretizations  are  motivated  by  the 
recent  theory  of  isogeometric  discrete  differential  forms  and  may  be  interpreted  as 
smooth  generalizations  of  Raviart-Thomas  elements.  The  new  discretizations  are  (at 
least)  patch-wise  C°  and  can  be  directly  utilized  in  the  Galerkin  solution  of  Darcy- 
Stokes-Brinkman  flow  for  single-patch  configurations.  When  applied  to  incompress¬ 
ible  flows,  these  discretizations  produce  pointwise  divergence-free  velocity  fields  and 
hence  exactly  satisfy  mass  conservation.  In  the  presence  of  no-slip  boundary  condi¬ 
tions  and  multi-patch  geometries,  the  discontinuous  Galerkin  framework  is  invoked 
to  enforce  tangential  continuity  without  upsetting  the  conservation  or  stability  prop¬ 
erties  of  the  method  across  patch  boundaries.  Furthermore,  as  no-slip  boundary 
conditions  are  enforced  weakly,  the  method  automatically  defaults  to  a  compatible 
discretization  of  Darcy  flow  in  the  limit  of  vanishing  viscosity.  The  proposed  dis¬ 
cretizations  are  extended  to  general  mapped  geometries  using  divergence-preserving 
transformations.  For  sufficiently  regular  single-patch  solutions,  we  prove  a  priori  er¬ 
ror  estimates  which  are  optimal  for  the  discrete  velocity  held  and  suboptimal,  by 
one  order,  for  the  discrete  pressure  held.  Our  estimates  are  additionally  robust  with 
respect  to  the  parameters  of  the  Darcy-Stokes-Brinkman  problem.  We  present  a  com¬ 
prehensive  suite  of  numerical  experiments  which  indicate  optimal  convergence  rates 
for  both  the  discrete  velocity  and  pressure  helds  for  general  configurations,  suggesting 
that  our  a  priori  estimates  may  be  conservative.  The  focus  of  the  current  paper  is 
strictly  on  incompressible  hows,  but  our  theoretical  results  naturally  extend  to  hows 
characterized  by  mass  sources  and  sinks. 

Keywords:  Darcy-Stokes-Brinkman  equations,  Generalized  Stokes  equations,  B-splines, 
Isogeometric  analysis,  Divergence-conforming  discretizations 
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1  Introduction 


The  Stokes  equations  describe  a  wide  variety  of  fluid  flows  where  advective  inertial 
forces  are  so  small  when  compared  with  viscous  forces  that  they  can  be  neglected 
altogether.  Such  flows  arise  in  a  large  number  of  applications  in  nature  and  technol¬ 
ogy,  from  the  flow  of  lava  in  the  Earth’s  mantle  [54]  to  microfluidic  flow  in  micro- 
electromechanical  devices  [41].  The  Darcy-Stokes-Brinkman  equations  are  a  simple 
extension  of  the  Stokes  equations  which  account  for  Darcy  drag  forces  in  highly  porous 
media  [13].  These  equations,  also  referred  to  as  the  generalized  Stokes  equations,  have 
been  used  to  model  groundwater  flow  [22],  heat  and  mass  transfer  in  pipes  [39],  and 
flow  in  biological  tissues  [40].  One  also  obtains  a  generalized  Stokes  problem  when 
nonlinear  terms  are  treated  explicitly  during  semi-implicit  time-integration  of  the 
unsteady  Navier-Stokes  equations. 

Despite  their  simple  appearance,  the  Stokes  and  generalized  Stokes  equations  have 
presented  considerable  difficulty  in  their  numerical  approximation.  At  the  heart  of 
the  matter  is  the  celebrated  Babuska-Brezzi  inf-sup  condition  [3,  12].  Simply  put,  the 
condition  states  that  one  must  properly  select  discrete  velocity  and  pressure  spaces 
in  order  to  arrive  at  a  stable  and  convergent  discrete  mixed  formulation.  Since  the 
inception  of  the  Marker-and-Cell  scheme  in  1965  [34],  a  large  number  of  finite  dif¬ 
ference,  finite  volume,  and  finite  element  methods  have  been  devised  to  address  the 
discrete  inf-sup  condition  in  the  context  of  the  Stokes  equations.  For  reference,  see 
the  review  by  Boffi,  Brezzi,  and  Fortin  [9] .  Most  methods  for  Stokes  flow  only  satisfy 
the  incompressibility  constraint  in  an  approximate  sense.  Some  bypass  the  inf-sup 
condition  entirely  through  the  use  of  a  stabilized  Petrov-Galerkin  method  [27].  How¬ 
ever,  methods  which  return  discretely  divergence-free  velocity  fields  are  generally  not 
robust  in  the  limit  of  vanishing  viscosity  when  applied  to  generalized  Stokes  flows 
[45].  Moreover,  mass  conservation  is  considered  to  be  of  prime  importance  for  cou¬ 
pled  flow-transport  [46],  and  it  has  been  demonstrated  that  methods  which  fail  to 
exactly  satisfy  the  incompressibility  constraint  suffer  from  spurious  velocity  oscilla¬ 
tions  when  applied  to  “high  pressure,  low  flow”  problems  [28,  44],  These  issues  have 
motivated  the  development  of  discretization  procedures  which  satisfy  the  incompress¬ 
ibility  constraint  exactly. 

One  of  the  simplest  methods  returning  a  divergence-free  velocity  field  is  the 
Pk  —  Pk~l  triangular  element  which  approximates  velocity  fields  using  continuous 
piecewise  polynomials  of  degree  k  and  pressure  fields  using  discontinuous  polynomi¬ 
als  of  degree  k  —  1.  This  method  satisfies  the  Babuska-Brezzi  condition  for  meshes 
containing  no  nearly-singular  vertices  provided  k  >  4  [53]  and  for  certain  macro¬ 
element  configurations  [2,  63].  Unfortunately,  the  method  is  not  stable  for  general 
meshes  and  polynomial  degrees.  Recently,  the  use  of  H(div;  f2)  elements  has  arisen 
as  a  new  paradigm  for  the  the  simulation  of  generalized  Stokes  flows  [38,  42,  61].  As 
these  approximations  are  typically  not  members  of  H1(D),  techniques  such  as  the 
interior  penalty  method  [1,  24,  62]  must  be  employed  to  enforce  tangential  continu- 
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ity  across  elements.  Alternatively,  one  can  modify  divergence-conforming  elements 
to  ensure  they  have  some  limited  tangential  continuity  [56],  and  some  authors  have 
elected  to  release  continuity  entirely  in  favor  of  hybridization  [18,  17].  In  the  same 
vein,  divergence-free  wavelets  have  been  proposed  for  the  solution  of  Stokes  flows 
[43,  58,  59,  55],  though  these  discretizations  have  a  complicated  construction  and  are 
entirely  limited  to  periodic  and  cuboidal  domains. 

In  this  paper,  we  present  new  divergence-conforming  B-spline  discretizations  for 
the  generalized  Stokes  problem.  These  discretizations  are  motivated  by  the  theory  of 
isogeometric  discrete  differential  forms  [15,  16]  and  may  be  interpreted  as  smooth  gen¬ 
eralizations  of  Raviart-Thomas  elements  [49].  The  new  discretizations  are  (at  least) 
patch-wise  C°  and  hence  can  be  directly  used  in  the  Galerkin  solution  of  generalized 
Stokes  flow  for  single-patch  configurations.  We  enforce  no-slip  boundary  conditions 
weakly  by  means  of  Nitsche’s  method  [48],  allowing  our  method  to  default  to  a  com¬ 
patible  discretization  of  Darcy  flow  in  the  limit  of  vanishing  viscosity.  Alternative 
inf-sup  stable  treatments  of  no-slip  boundary  conditions  have  been  investigated  in 
the  context  of  Stokes  flow  in  [14].  In  the  presence  of  multi-patch  geometries,  we  in¬ 
voke  the  discontinuous  Galerkin  framework  in  order  to  enforce  tangential  continuity 
across  patch  interfaces  while  maintaining  the  stability  and  conservation  properties 
of  the  method.  The  proposed  discretizations  are  extended  to  general  mapped  ge¬ 
ometries  using  divergence-  and  integral-preserving  transformations.  For  single-patch 
solutions,  we  are  able  to  prove  robust  a  priori  error  estimates  which  are  optimal 
for  the  discrete  velocity  held  and  suboptimal,  by  one  order,  for  the  discrete  pressure 
held.  This  is  reminiscent  of  error  estimates  for  stabilized  equal-order  interpolations 
of  the  Stokes  equations  [27,  37].  Our  derived  estimates  are  also  robust  with  respect 
to  the  parameters  of  generalized  Stokes  how.  We  utilize  the  methods  of  exact  and 
manufactured  solutions  to  validate  our  estimates,  and  we  find  that  pressure  actually 
converges  at  optimal  order.  We  further  test  the  effectiveness  of  our  method  using  a 
collection  of  benchmark  problems:  two-dimensional  creeping  lid-driven  cavity  how, 
three-dimensional  creeping  lid-driven  cavity  how,  and  Darcy-dominated  generalized 
Stokes  how  subject  to  boundary  layers. 

An  outline  of  this  paper  is  as  follows.  In  the  following  section,  we  present  some 
basic  notation.  In  Section  3,  we  recall  the  generalized  Stokes  problem  subject  to  ho¬ 
mogeneous  boundary  conditions.  In  Section  4,  we  briefly  review  B-splines,  the  basic 
building  blocks  of  our  new  discretization  technique,  and  in  Section  5,  we  define  the  B- 
spline  spaces  which  we  will  utilize  to  discretize  velocity  and  pressure  helds.  In  Section 
6,  we  present  our  discrete  variational  formulation  for  the  generalized  Stokes  problem 
and  prove  continuity,  stability,  and  a  priori  error  estimates  for  the  single-patch  set¬ 
ting.  In  Section  7,  we  discuss  how  to  extend  our  methodology  to  the  multi-patch 
setting.  In  Sections  8  and  9,  we  present  numerical  results,  and  in  Section  10,  we  draw 
conclusions.  Throughout  this  paper,  we  make  explicit  all  of  our  estimates’  dependen¬ 
cies  on  the  problem  parameters  as  well  as  the  penalty  parameter  of  Nistche’s  method. 
This  will  require  a  lengthy  and  tedious  analysis,  but  we  believe  that  knowledge  of 
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such  dependencies  is  of  great  practical  importance.  Additionally,  the  focus  of  this 
paper  is  strictly  on  incompressible  flows.  Our  theoretical  results  naturally  extend  to 
flows  characterized  by  mass  sources  and  sinks. 

2  Notation 

We  begin  this  paper  with  some  basic  notation.  For  d  a  positive  integer  representing 
dimension,  let  D  C  denote  an  arbitrary  bounded  Lipschitz  domain  with  boundary 
d D.  As  usual,  let  L2(D)  denote  the  space  of  square  integrable  functions  on  D  and 
define  L 2(D)  =  ( L2(D))d .  We  will  also  utilize  the  more  general  Lebesgue  spaces 
LP(D)  where  1  <  p  <  oo  and  their  vectorial  counterpart  ~LP(D).  Let  Hk\D )  denote 
the  space  of  functions  in  L2(D )  whose  kth- order  derivatives  belong  to  L2(D)  and 
define  Hk(D)  =  ( Hk(D))d .  We  identify  with  Hk(D )  the  standard  Sobolev  norm 

\\V\\ H*(D)  =  (  H-D<*VH L2(D)  ) 

\M<*  / 

where  a  =  (a±,  a2,  ■  ■  ■ ,  eta)  is  a  multi-index  of  non- negative  integers,  |a|  =  ot\  +  a2  + 
. . .  +  aa,  and 

Q\<*\ 

D  =  dxTdx*2 . . .  dxQdd ' 

We  denote  the  Sobolev  semi-norms  as  |  •  \h*(d)i  and  we  use  the  convention  H°(D)  = 
L2(D).  Throughout,  Sobolev  spaces  of  fractional  order  are  defined  using  function 
space  interpolation  (see,  e.g.,  Chapter  1  of  [57]).  We  define  H(\  (D)  c  Hl(D)  to  be 
the  subspace  of  functions  with  homogeneous  boundary  conditions  and  define  Hq  (D) 
to  be  the  vectorial  counterpart  of  Hq(D).  We  define  Hs(div;D)  to  be  the  Sobolev 
space  of  all  functions  in  Hs(fl)  whose  divergence  also  belongs  to  HS(D).  This  space 
is  equipped  with  the  norm 

1  /9 

||v||H*(div;r>)  =  (IIvIIhs(d)  +  lldivvlllqD))  • 

When  s  =  0,  we  drop  the  index.  We  also  define 

H0(div;  D)  =  {v  e  H(div;  D)  :  v  •  n  =  0  on  d D} 

where  n  denotes  the  outward  pointing  unit  normal.  Finally,  we  denote  Lq(D)  C  L2(D) 
as  the  space  of  square- integrable  functions  with  zero  average  on  D. 

3  The  Generalized  Stokes  Problem 

In  this  section,  we  recall  the  generalized  Stokes  problem  subject  to  homogeneous 
Dirichlet  boundary  conditions.  For  d  a  positive  integer,  let  denote  a  Lipschitz 
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bounded  open  set  of  Md.  Throughout  this  paper,  d  will  be  either  2  or  3.  The  problem 
of  interest  is  as  follows. 

(  Given  a  :  12  — >  R,  v  :  12  — >  R,  and  f :  12  — >  Rd,  find  u  :  12  — >■  Rd  and  p  :  12  — >  R 
such  that 


(S)< 


tru  —  V  ■  (2z/Vsu)  +  gradp  = 

f 

in 

(1) 

divu  = 

0 

in  12 

(2) 

u  = 

0 

on  (912. 

(3) 

Above,  u  denotes  the  flow  velocity  of  a  fluid  moving  through  the  domain  12,  p  denotes 
the  pressure  acting  on  the  fluid  divided  by  the  fluid  density,  v  denotes  the  kinematic 
viscosity  of  the  fluid,  o  denotes  the  reaction  coefficient  which  gives  the  ratio  of  the 
viscosity  to  the  permeability  of  the  fluid,  and  f  denotes  a  body  force  acting  on  the 
fluid  divided  by  the  density.  We  assume  that  the  viscosity  is  taken  to  be  uniformly 
positive  (i.e.,  3z/0  >  0  such  that  v  >  z/0)  and  that  the  reaction  coefficient  is  taken  to 
be  non- negative  (i.e.,  a  >  0).  Note  that  the  pressure  is  only  determined  up  to  an 
arbitrary  constant. 

Let  us  now  make  an  assumption  regarding  the  data  of  our  problem.  Notably,  let 
us  assume  that  cr,  v  G  L°°(12)  and  f  G  L3(12).  The  weak  form  for  the  generalized 
Stokes  problem  is  then  written  as  follows: 


(  Find  u  G  and  p  G  such  that 


a(u,  v)  -  b(p,  v)  +  b(q,  u)  =  (f,  v)La(n) 


(4) 


(W){ 


for  all  v  G  and  q  G  Lg(fl)  where 


a(w,v)  =  (2z/Vsw,  Vsv)(L2(n))dxd  +  (aw,v)L2(n)  ,  Vw,vgHJ(0),  (5) 
%,v)  =  (q,  divv)i2(a),  Vq  G  Lg(fl),v  G  Ho(fl).  (6) 


We  have  the  following  theorem  which  is  a  simple  consequence  of  coercivity,  continuity, 
and  a  continuous  inf-sup  condition. 

Theorem  3.1.  Problem  ( W )  has  a  unique  weak  solution  (u ,p)  G  Hg(f2)  x«. 


4  B-splines  and  Geometrical  Mappings 

In  this  section,  we  briefly  introduce  B-splincs,  the  primary  ingredient  in  our  dis¬ 
cretization  technique  for  the  generalized  Stokes  equations.  We  also  introduce  map¬ 
pings  which  will  allow  us  to  extend  our  discretization  technique  to  general  geometries 


5 


of  engineering  interest.  For  an  overview  of  B-splines,  their  properties,  and  robust  al¬ 
gorithms  for  evaluating  their  values  and  derivatives,  see  de  Boor  [23]  and  Schumaker 
[52],  For  the  application  of  B-splines  to  finite  element  analysis,  see  Hollig  [35]  and 
Cottrell,  Hughes,  and  Bazilcvs  [20]. 

4.1  Univariate  B-splines 

For  two  positive  integers  k  and  n,  representing  degree  and  dimensionality  respectively, 
let  us  introduce  the  ordered  knot  vector 

S:={0  =  £i, (7) 


where 


Cl  —  £2  5;  •  •  •  Cn+fc+1- 


Given  H  and  k ,  univariate  B-splinc  basis  functions  are  constructed  recursively  starting 
with  piecewise  constants  (k  —  0): 


B°(0 


1  if  6  <  C  <  6+1 

0  otherwise. 


(8) 


For  k  —  1,  2,  3, . . .,  they  are  defined  by 


= 


C~6 

6+fc  6 


+ 


6 


i+k+ 1 


C 


6+fc+i  6+1 


(9) 


When  6+ k  —  6  —  0,  ^  %  is  taken  to  be  zero,  and  similarly,  when  6+fc+i  —  6+1  —  0, 

S  i+k  S  i 

$i+fc+i-g  js  taken  to  be  zero.  B-spline  basis  functions  are  piecewise  polynomials  of 

s  i  “F  ^  “I- 1 

degree  k,  form  a  partition  of  unity,  have  local  support,  and  are  non-negative.  An 
example  of  a  cubic  B-spline  basis  is  shown  in  Fig.  1.  Note  the  basis  is  C 2  everywhere 
in  the  interval  (0,1).  Enhanced  smoothness  is  one  of  the  defining  features  of  B- 
splines.  We  refer  to  linear  combinations  of  B-spline  basis  functions  as  B-splines  or 
simply  splines. 

Let  us  now  introduce  the  vector  £  =  {£1, . . . ,  Cm}  of  knots  without  repetitions  and 
a  corresponding  vector  {ri, . . . ,  rm}  of  knot  multiplicities.  That  is,  rt  is  defined  to  be 
the  multiplicity  of  the  knot  Q  in  5.  By  construction,  Y77=\  ri  —  n  +  k  + 1.  We  assume 
that  rl  <  k+1.  Let  us  further  assume  throughout  that  ry  =  rm  =  k+ 1,  i.e,  that  5  is  an 
open  knot  vector.  This  allows  us  to  easily  prescribe  Dirichlet  boundary  conditions. 
At  the  point  Q,  B-spline  basis  functions  have  ay  k  —  ry  continuous  derivatives. 
Therefore,  —  1  <  ay  <  p—  1,  and  the  maximum  multiplicity  allowed,  ry  =  k  +  1,  gives 
a  discontinuity  at  Q.  We  define  the  regularity  vector  a  by  a  :=  (ay, . . .  ,am}.  By 
construction,  a  1  =  am  =  —1.  In  what  follows,  we  utilize  the  notation 


a  |  =  minjay  :  2  <  i  <  m  —  1} 


(10) 
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Figure  1:  Cubic  B-spline  basis  functions  for  open,  non-uniform  knot  vector  5  = 
{0,  0,  0,  0, 1/3, 2/3, 1, 1, 1, 1}.  Note  the  basis  is  C2  everywhere  in  the  interval  (0, 1). 

and  a  —  1  :=  {—1,  ol-i  —  1, . . . ,  am_i  —  1,  —1}  when  ai  >  0  for  2  <  i  <  m  —  1. 

We  denote  the  space  of  B-splines  spanned  by  the  basis  functions  £>/  as 

Sy=span{B‘}’l1.  (11) 

When  k  >  1  and  CKj  >  0  for  2  <  i  <  m  —  1,  the  derivatives  of  functions  in  S are 
splines  as  well.  In  fact,  we  have  the  stronger  relationship 

<i2> 

That  is,  the  derivative  operator  dx  :  — )■  S'^T1 1  is  surjective.  One  of  the  most 

important  properties  of  univariate  B-splines  is  refinement  and,  perhaps  more  impor¬ 
tantly,  nestedness  of  refinement.  Notably,  knot  insertion  and  degree  elevation  allow 
one  to  define  a  sequence  of  nested,  refined  B-spline  spaces.  Knot  insertion  and  degree 
elevation  algorithms  are  described  in  detail  in  Chapter  2  of  [20]. 

4.2  Multivariate  B-splines 

The  definition  of  multivariate  B-splines  follows  easily  through  a  tensor-product  con¬ 
struction.  For  d  again  a  positive  integer,  let  us  consider  the  unit  cube  12  =  (0,  l)d  C 
Md,  which  we  will  henceforth  refer  to  as  the  parametric  domain.  Mimicking  the  one¬ 
dimensional  case,  given  integers  ki  and  n/  for  Z  =  1,  ...,d,  let  us  introduce  open 
knot  vectors  Et  =  {£M, . . . ,  t,ni+kl+i,i}  and  the  associated  vectors  Ci  =  {Ci ,h  ■  ■  ■ ,  C mlti}, 
{ri  ;, . . . ,  rmhi},  and  oli  =  {aq,;, . . . ,  <ymhi}.  There  is  a  parametric  Cartesian  mesh  M.h 


7 


associated  with  these  knot  vectors  partitioning  the  parametric  domain  into  rectangu¬ 
lar  parallelepipeds.  Visually, 


■M-h  —  {Q  —  (C iuh  C ii+l,l)  j  1  <  k  <  mi  ~  1}  • 


(13) 


For  each  element  Q  G  A4h  we  associate  a  parametric  mesh  size  Jiq  =  diarn(Q).  We 
also  dehne  a  shape  regularity  constant  A  which  satisfies  the  inequality 

A-i  <  <  A,  VQ  e  (14) 

llQ 

where  hq^ in  denotes  the  length  of  the  smallest  edge  of  Q.  A  sequence  of  parametric 
meshes  that  satisfy  the  above  inequality  for  an  identical  shape  regularity  constant  is 
said  to  be  locally  quasi-uniform. 

We  associate  with  each  knot  vector  S;  (/  =  1 , ,d)  univariate  B-spline  basis 
functions  B^\  of  degree  /q  for  h  =  1, . . . ,  rq.  On  the  mesh  M.h,  we  dehne  the  tensor- 
product  B-spline  basis  functions  as 

Biu---fdd  :=  Bh,  i  ®  •  •  ■  0  Bild >  *i  =  ...  (15) 


We  then  accordingly  dehne  the  tensor-product  B-spline  space  as 


=  S*\';;;:%d(Mh)  :=  span 


. . na 

” . J  i,.i . «=i 


(16) 


The  space  is  fully  characterized  by  the  mesh  At/,,,  the  degrees  ki ,  and  the  regularity 
vectors  ex;,  as  the  notation  rehects.  Like  their  univariate  counterparts,  multivariate 
B-spline  basis  functions  are  piecewise  polynomial,  form  a  partition  of  unity,  have  local 
support,  and  are  non-negative.  Defining  the  regularity  constant 


a  \=  min  min  (17) 

l=l,...,d2<ii<mi-l 

we  see  that  our  B-splines  are  (^"-continuous  throughout  the  domain  Q.  Refinement  of 
multivariate  B-spline  bases  is  obtained  by  applying  knot  insertion  and  degree  elevation 
in  tensor-product  fashion.  In  the  remainder  of  the  text,  we  consider  a  family  of 
nested  meshes  {. Mh}h<ho  and  associated  B-spline  spaces  h<h  that 

have  been  obtained  by  successive  applications  of  knot  refinement.  Furthermore,  we 
assume  throughout  that  the  mesh  family  {Af/i}fe</lo  is  locally  quasi-uniform. 

Note  that  each  element  Q  =  <S>i=i,...,d  (C iuh  Cn+i,i)  das  the  equivalent  representation 
Q  =  ®i=i,...,d  Cji+i,i)  f°r  some  index  ji.  With  this  in  mind,  we  associate  with  each 
element  a  support  extension  Q,  defined  as 

Q  :=  ®Z=l,...,d  £ji+Pi+l,l)  ■  (18) 

The  support  extension  is  the  interior  of  the  set  formed  by  the  union  of  the  supports 
of  all  B-spline  basis  functions  whose  support  intersects  Q.  Note  that  each  element 
belongs  to  the  support  extension  of  at  most  Ih=lv  jrf(2pz  +  1)  elements.  The  support 
extension  is  a  natural  object  to  consider  when  examining  the  local  approximation 
properties  of  a  B-spline  space. 


4.3  Piecewise  Smooth  Functions,  Geometrical  Mappings,  and 
Physical  Mesh  Entities 

On  the  parametric  mesh  Aih,  we  define  the  space  of  piecewise  smooth  functions  with 
interelement  regularity  given  by  the  vectors  op, . . . ,  as 


/^OO 

^CKi, 


_  OO 

,OLd  ^(Xl,...,OLd 


(. M , 


(19) 


Precisely,  a  function  in  C™  ad  is  a  function  whose  restriction  to  an  element  Q  G  Ain 
admits  a  C°°  extension  in  the  closure  of  that  element  and  which  has  atlj  contin¬ 
uous  derivatives  with  respect  to  the  Zth  coordinate  along  the  internal  mesh  faces 
{(xi,...,a;d)  :  xt  =  Quh  QVtV  <  xv  <  ^  0  for  all  ix  =  2,...,mz  -  1  and 

ji /  =  1, . . .  ,mi>  —  1.  Note  immediately  that  any  function  lying  in  the  B-spline  space 
also  lies  in 

Unless  specified  otherwise,  we  assume  throughout  the  rest  of  the  paper  that  the 
physical  domain  fi  C  I6*  can  be  exactly  parametrized  by  a  geometrical  mapping 
F  :  D  — >  D  belonging  to  (C^j  )  with  piecewise  smooth  inverse.  We  further 
assume  that  the  physical  domain  is  simply  connected  with  connected  boundary 
dfl  and  the  geometrical  mapping  is  independent  of  the  mesh  family  index  h.  See,  for 
example,  the  sequence  of  mapped  meshes  depicted  in  Figure  2.  A  geometrical  mapping 
meeting  our  criteria  could  be  defined  utilizing  B-splines  or  Non-Uniform  Rational  B- 
Splines  (NURBS)  on  the  coarsest  mesh  A4/,0.  For  examples  of  such  mappings,  see 
Chapter  2  of  [20].  NURBS  mappings  are  especially  useful  as  they  can  represent  many 
geometries  of  scientific  and  engineering  interest  and  are  the  main  tools  employed 
in  Computer  Aided  Design  (CAD)  software.  Later  in  this  paper,  we  will  utilize 
polar  mappings  to  solve  flow  problems  on  cylindrical  geometries  in  order  to  preserve 
symmetries. 

The  geometrical  mapping  F  naturally  induces  a  mesh 


Kh  =  {K  :K  =F(Q),QeMh} 


(20) 


on  the  physical  domain  D.  We  define  for  each  element  K  G  a  physical  mesh  size 


liR  —  ll-DF \\L°°{Q)hQ  (21) 

where  Q  is  the  pre-image  of  K ,  and  we  also  define  the  support  extension  K  =  F(Q). 
We  define  for  a  given  mesh  the  global  mesh  size 


h  =  max  {hx,  K  G  /C^}  . 

Note  that  as  the  parametric  mesh  family  {A ^h}h<h0  locally  quasi-uniform  and  the 
geometrical  mapping  F  is  independent  of  the  mesh  family  index  h,  the  physical  mesh 
family  {^h}h<h0  afo°  l°caUy  quasi-uniform.  We  refer  to  the  physical  domain  D  and 
its  pre-image  D  interchangeably  as  the  patch.  It  should  be  noted  that,  in  general,  the 
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Coarse  Mesh 


First  Refinement 


Second  Refinement 

Figure  2:  Illustration  in  the  two-dimensional  setting  of  how  the  parametric  mapping 
F  is  independent  of  the  mesh  family  index  h. 
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domain  fl  cannot  be  represented  using  just  a  single  patch.  Instead,  multiple  patches 
must  be  employed.  We  will  discuss  further  the  multi-patch  setting  in  Section  7. 

We  define  on  the  parametric  mesh  a  set  of  mesh  faces  Fh  =  {F}  where  F  is  a  face 
of  one  or  more  elements  in  A ip-  We  define  the  physical  set  of  mesh  faces  as 

Fh  =  {F=F(F):FeFh} 
and  we  define  the  boundary  mesh  to  be 

Th  =  {F  e  Fh  :  F  C  ,911}  . 


By  construction, 


<911  —  Up(zrhF- 


Note  that  for  each  face  F  G  Th  there  is  a  unique  K  e  JCh  such  that  F  is  a  “face”  of 
K  (in  the  sense  that  F  is  the  image  of  a  face  of  Q,  the  pre-image  of  K).  We  hence 
define  for  such  a  face  the  mesh  size 


hp  :=  hjp . 

One  may  also  define  hp  to  be  the  wall-normal  mesh-size  as  is  done  in  [7].  Such  a 
definition  is  more  appropriate  when  stretched  meshes  are  utilized  in  the  presence  of 
layers. 

Throughout  the  paper,  we  will  utilize  the  terminology  “a  constant  independent 
of  IF.  When  we  employ  such  terminology,  we  simply  indicate  that  the  constant 
will  not  depend  on  the  given  mesh  and,  in  particular,  its  size.  The  constant  may, 
however,  depend  on  the  domain,  the  shape  regularity  of  the  parametric  mesh  family, 
the  polynomial  degrees  of  the  employed  B-splinc  spaces,  and  global,  mesh-invariant 
measures  of  the  parametric  mapping. 


5  Discretization  of  Velocity  and  Pressure  Fields 

In  this  section,  we  define  the  B-spline  spaces  which  we  will  utilize  to  discretize  the 
velocity  and  pressure  fields  appearing  in  the  generalized  Stokes  problem.  These  spaces 
are  motivated  by  the  recent  theory  of  isogeometric  discrete  differential  forms  [15,  16] 
and  may  be  interpreted  as  smooth  generalizations  of  Raviart-Thomas  elements  [49] . 
We  first  define  our  discrete  velocity  and  pressure  spaces  on  the  parametric  domain 
fl  =  (0,  l)d  and  then  define  discrete  spaces  on  the  physical  domain  using  divergence- 
and  integral-preserving  transformations.  We  finish  this  section  with  a  presentation 
of  local  approximation  estimates  and  trace  inequalities  for  our  discrete  velocity  and 
pressure  spaces. 
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5.1  Discrete  Spaces  on  the  Parametric  Domain 

Using  the  notation  of  the  previous  section  and  assuming  that 

a  :=  min{|a;|  :  l  —  1, . . . ,  d}  >  1, 


we  define  the  following  two  spaces: 


Vh:= 


oki,k2~l  CKi~l,K2 

^  Otl ,0:2  —  1  A  ^CKi  — 1,0-2 

Qki  —  I,k2,k3  —  1  Qki  —  l,k2  —  l,k3 


'ik±  —  l,k2 


qki,k2  —  l,k3  —  1  w  Q/Ci  —  1,/C2,«3  —  1  v  Q/e  1  — 1,AC2  — 1,^3  j  q 

^011,0:2  —  1,0:3  —  1  *  ^01  —  1,02,03—  1  ^01—1,02  —1,03  11  11'  D, 


if  d  —  2, 


2-) J  °c*i-l,c*2-l  11  Ll  —  A 

-^-h  ■  \  qki  —  l,fc2  —  l,fc3  —  1  ;r  7  o 

L  °CK1-1, 012  —  1,  «3  — 1  11  a  — 

The  space  Vh  comprises  our  set  of  discrete  velocity  fields  while  Qh  comprises  our 
set  of  discrete  pressure  fields.  Note  that  as  cc  >  1,  our  discrete  velocity  fields  are 
H1-conforming.  If  we  allow  a  =  0,  our  spaces  collapse  to  standard  Raviart- Thomas 
mixed  finite  elements  [49].  In  order  to  deal  with  no-penetration  boundary  conditions, 
we  make  use  of  the  following  constrained  discrete  spaces: 

Vo ,h  ■=  { G  Vh  :  v,,.  •  n  =  0  on  dfl)  , 


Qo,h  ■=  e  Qh  :  J  qhdx  =  0  j  . 

Above,  n  denotes  the  outward-facing  normal  to  <9fh  As  specified  in  the  introduction, 
we  choose  to  enforce  no-slip  boundary  condition  weakly  using  Nitsche’s  method  [48]. 
Due  to  the  special  relationship  given  by  (12),  the  spaces  Vo x  and  <20,h  along  with  the 
parametric  divergence  operator  form  the  bounded  discrete  cochain  complex 


where  div  is  the  divergence  operator  on  the  unit  cube  fh  In  fact,  we  have  a  much 
stronger  result  due  to  the  results  of  [15]. 

Proposition  5.1.  There  exist  L2- stable  projection  operators  flk  :  Ho(div;fl)  — )■  Vo ,h 
and  flT  :  L'l(Q)  — >  Qo  h  such  that  the  following  diagram  commutes: 

55 lh 


H0(div;  Q) 


Furthermore,  there  exists  a  positive  constant  Cu  independent  of  h  such  that 
II  v||Hi(n)  <  Cu||v||Hl(fi),  Vv  G  H0(div;  Q)  n  H1^). 
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5.2  Discrete  Spaces  on  the  Physical  Domain 

To  define  our  discrete  velocity  and  pressure  spaces  on  the  physical  domain,  we  intro¬ 
duce  the  following  pullback  operators: 

tu(v)  :=  det(DF)(DF)_1(voF),  vGH0(div;D)  (24) 

iP{q)  :=  det  (DF)  (q  o  F) ,  q  G  Lq(Q)  (25) 

where  -DF  is  the  Jacobian  matrix  of  the  parametric  mapping  F.  The  push-forward 
given  by  (24),  popularly  known  as  the  Piola  transform,  has  two  important  properties: 
(i)  it  preserves  the  nullity  of  normal  components,  (ii)  it  maps  divergences  to  diver¬ 
gences.  Hence,  it  maps  divergence-free  fields  in  parametric  space  to  divergence-free 
fields  in  physical  space,  as  illustrated  in  Fig.  3.  On  the  other  hand,  the  push-forward 
given  by  (25)  has  the  property  that  it  preserves  the  nullity  of  the  integral  operator. 
Due  to  these  properties,  we  have  the  following  commuting  diagram: 


H0(div;O)  -^4  Lg(fi) 


lu  ip  (26) 

H0(div;O)  -^4  Lg(fi). 

This  motivates  the  use  of  the  following  discrete  velocity  and  pressure  spaces  in  the 
physical  domain: 

V0 ,h  :=  jv  6  H0(div;  Q)  :  iu(v)  G  V0)/l|  , 


Qo,h  {<?  £  I/q(D)  :  tp(q)  G  Qo,h|  • 

Furthermore,  we  define  the  projectors  :  H0(div;O)  — >  Vo,h  and  n qh  :  L2(Q)  — > 
Qo,h  via  the  compositions 


nL 


:=  t.. 1  ° 


n°4 


O  K 


UQh 


:  =  1  o 


n°4 


O  t 


v 


Employing  the  preceding  results  and  terminology  as  well  as  the  smoothness  properties 
of  the  parametric  mapping  F ,  we  arrive  at  the  following  proposition. 


Proposition  5.2.  The  following  diagram  commutes: 


H0(div;  Q) 


div 


»  Lim 


n° 


ft0- 

Qh 


(27) 


V, 


div 


0  ,h 


t  Qo,h- 


Furthermore,  there  exists  a  positive  constant  Cu  independent  of  h  such  that 


in0 


<  Cu||  v|| Hi(Q} ,  Vv  e  H0(div;  Q)  fl  H1(0). 


(28) 
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Figure  3:  The  Piola  transform  maps  divergence-free  fields  in  parametric  space  to 
divergence-free  fields  in  physical  space,  as  shown  here  for  the  case  of  a  divergence-free 
B-spline. 


We  immediately  have  an  inf-sup  condition  for  onr  discrete  velocity /pressure  pair. 

Proposition  5.3.  There  exists  a  positive  constant  f3  independent  of  h  such  that  the 
following  holds:  for  every  qh  E  Qo,h,  there  exists  a  E  Vo,h  such  that: 

divv/j  =  qh  (29) 


and 


Hence, 


Iv/iIIh1^)  <  $  lk/i||z,2(n). 


(divv/i,  qh)L2{ti) 

ih&Qo,h  VhGv0,h  llvftllHhmlMU2^) 
qh¥=  o 


inf 


sup 


>»■ 


(30) 

(31) 


Proof.  Let  qh  G  Qo,h  be  arbitrary.  It  is  a  classical  result  (see  pg.  24  of  [29],  for 
example)  that  there  exists  a  function  v  G  Hq(O)  such  that  divv  =  qh  and 

IMIh1^)  <  ^_1|I^I|l2(0) 

where  (3  is  a  positive  constant  independent  of  v.  Let  =  Lly  v.  Then,  by  Proposi¬ 
tion  5.2,  divv/,  =  divlly  v  =  Llgh  divv  =  qh  and 

1 1  1 1 H1  (ri)  —  ^ti||  v||  H1(f2)  —  Cuf3  1||0,h|U2(Q)- 

Thus  the  theorem  holds  with  B  —  ff.  □ 
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We  would  like  to  state  that  such  an  inf-sup  condition  would  not  have  held  if  we 
had  elected  to  strongly  enforce  homogeneous  tangential  boundary  conditions  within 
the  space  Vo,/,.  This  is  because 

div{v0,fenHj(n)}  c  Q0)h. 

Hence,  there  exists  a^e  Qo,/i  such  that  qh  0  and 

gu  (divvfe,  gh)L2(n)  = 

v^eVo^nH^Q)  II II H1(f2) hhWmn) 

Hence,  an  alternative  methodology  must  be  employed  if  one  wishes  to  strongly  enforce 
homogeneous  tangential  boundary  conditions  within  the  space  Vo,/,.  For  example,  in 
[14],  a  special  discrete  pressure  space  was  constructed  in  the  two-dimensional  setting 
by  selectively  reducing  the  dimensionality  of  Qo  h  using  T-splines  [5].  However,  a  proof 
of  mesh-independent  discrete  stability  remains  absent  with  this  choice  of  pressure 
space,  and  the  convenient  tensor-product  structure  of  B-splines  is  lost. 

We  also  have  the  following  result. 

Proposition  5.4.  7/v/,  G  Vo,/,  satisfies 

(di Wft,9ft)L2([J)  =  0,  Wqh  G  Q0,h,  (32) 


then  divv/,  =  0. 


Proof.  The  proof  holds  trivially  as  div  maps  Vo./,  onto  □ 

Hence,  by  choosing  Vo,/,  and  Qo,h  as  discrete  velocity  and  pressure  spaces,  we 
arrive  at  a  discretization  that  automatically  returns  velocity  fields  that  are  pointwise 
divergence-free. 

5.3  Visualization  of  Basis  Functions 

We  now  visualize  some  of  the  basis  functions  associated  with  our  chosen  discrete 
velocity  and  pressure  spaces.  For  ease  of  presentation,  we  confine  ourselves  to  the 
two-dimensional  setting.  We  also  work  with  the  unconstrained  analogues  of  Vo,/,  and 
Qo,h,  which  we  denote  as  V/,  and  Q/,.  Let  k\  —  k2  —  2,  and  let  Eq  and  S2  be  equal  to 

Si  :=S2  :=  {0,0,  0,1/3,  2/3, 1,1,1}. 

These  polynomial  degrees  and  knot  vectors  define  a  parametric  mesh  At/,  and  B-spline 
spaces  V/,  and  <2/,  over  A4h.  To  define  the  physical  domain,  we  employ  a  biquadratic  B- 
spline  mapping.  The  control  net  defining  this  mapping  (see  Chapter  2  of  [20])  and  the 
resulting  physical  mesh  /C/,  are  illustrated  in  Figure  4.  In  Figure  5,  we  have  depicted 
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a  “first-component”  vector  basis  function  of  the  discrete  velocity  space  1 4-  Note 
that  the  basis  function  is  (^-continuous  along  horizontal  parametric  lines  and  C°- 
continuous  along  vertical  parametric  lines.  Further  note  that  the  directionality  of  the 
basis  function  is  preserved  under  the  map  iu  in  the  sense  that  the  function  is  oriented 
in  the  direction  of  horizontal  parametric  lines  in  both  parametric  and  physical  space. 
In  Figure  6,  we  have  depicted  a  typical  “second-component”  vector  basis  function  of 
14.  Note  that  the  basis  function  is  C°-continuous  along  horizontal  parametric  lines 
and  (^-continuous  along  vertical  parametric  lines,  and  the  directionality  of  the  basis 
function  is  preserved  under  the  map  iu  in  the  sense  that  the  function  is  oriented  in  the 
direction  of  vertical  parametric  lines  in  both  parametric  and  physical  space.  Finally, 
in  Figure  7,  we  have  depicted  a  representative  basis  function  of  the  discrete  pressure 
space  Qh  which  is  C°-continuous  in  both  parametric  and  physical  space. 

5.4  Approximation  Results  and  Trace  Inequalities 

Let  us  define 

k' —  min  \pi  —  1| .  (33) 

Note  that  the  discrete  velocity  and  pressure  spaces  Voy  and  Q0,/)  consist  of  mapped 
piecewise  polynomials  which  are  complete  up  to  degree  k' .  Hence,  k'  may  be  thought 
of  as  the  polynomial  degree  of  our  discretization  technique.  The  following  result 
details  the  local  approximation  properties  of  our  discrete  spaces.  Its  proof  may  be 
found  in  [15]. 

Proposition  5.5.  Let  K  e  /Q,  and  K  denote  the  support  extension  of  K .  For 
0  <  j  <  s  <  k'  +  1,  we  have 

|v-n^v|Hj(K)<C/>y-*||v||H.(#),  VveH*(A')nH„(div;f!)  (34) 
k  -  <  Ch^>\\q\\H,{k),  v9  e  H’(k)  n  ig(fi)  (35) 

where  C  denotes  a  positive  constant,  possibly  different  at  each  appearance,  independent 
ofh. 

Hence,  our  discrete  spaces  deliver  optimal  rates  of  convergence  from  an  approx¬ 
imation  point  of  view.  This  being  said,  the  results  of  Proposition  5.5  are  riddled 
with  inconvenient  interpolation  constants  C  which  depend  on,  among  other  things, 
the  polynomial  degree  and  continuity  of  our  approximation  spaces.  To  attack  the 
question  of  degree  and  continuity  directly,  Bcirao  da  Veiga  et  al.  derived  interpola¬ 
tion  estimates  for  B-splines  with  explicit  dependence  on  degree  and  continuity  in  [21]. 
However,  the  derived  estimates  are  only  available  for  interpolations  of  Hcrmite-type. 
All  of  this  seems  to  indicate  that  function  analytic  estimates  have  their  limitations. 
Alternatively,  one  can  use  numerics  to  study  the  approximation  properties  of  discrete 
spaces  using  the  theory  of  Kolmogorov  n-widths.  This  approach  allows  one  to  exactly 
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Control  Mesh 


Physical  Mesh 

Figure  4:  The  control  net  and  physical  mesh  for  the  biquadratic  B-spline  surface  with 
S,  :=S2:=  {0,0,  0,1/3,  2/3, 1,1,1}. 
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Figure  5:  Vector  plots  of  a  representative  first-component  vector  basis  function  of  the 
discrete  velocity  space  Vh  in  both  parametric  and  physical  space. 
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Figure  6:  Vector  plots  of  a  representative  second-component  vector  basis  function  of 
the  discrete  velocity  space  Vh  in  both  parametric  and  physical  space. 
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Parametric  Space 


Physical  Space 

Figure  7:  Contour  plots  of  a  representative  basis  function  of  the  discrete  pressure 
space  Qh  in  both  parametric  and  physical  space. 
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compute  the  interpolation  constants  associated  with  variational  projection  through 
the  solution  of  a  generalized  eigenproblcm.  The  theory  of  Kolmogorov  n-widths  was 
used  to  study  one- dimensional  B-spline  discretizations  in  [25].  This  paper  revealed 
that  maximal  continuity  B-spline  spaces  harbor  nearly  optimal  resolution  properties 
and  admit  smaller  intepolation  constants  than  lower  continuity  spaces.  Recently, 
the  theory  of  Kolmogorov  n-widths  has  been  used  to  study  multi-dimensional  and 
compatible  B-spline  discretizations.  This  study  has  also  revealed  the  advantage  of 
employing  B-splines  of  maximal  continuity  in  the  multi- dimensional  setting.  The 
results  of  this  study  will  be  presented  in  a  forthcoming  publication. 

We  will  also  need  the  following  trace  estimate  in  our  proceeding  mathematical 
analysis.  Its  proof  can  be  found  in  [26]. 

Proposition  5.6.  Let  K  G  /C*  and  Q  =  F_1(K').  Then  we  have 

||  (Vsv/j)  n||(i2 {dK))d  <  Ctracehf} 1 1 v/j 1 1 Hi (x-} ,  VvA  G  V0, h  (36) 

where  C trace  denotes  a  positive  constant  independent  of  h. 

In  [26],  it  was  shown  that  Proposition  5.6  holds  with  Ctrace  ~  (k')2.  However,  our 
numerical  experience  has  indicated  that  a  corresponding  global  trace  inequality  holds 
With  Ctrace  n-u  k'  if  B-splines  of  maximal  continuity  are  utilized.  This  allows  us  to 
select  a  smaller  penalty  parameter  when  employing  Nitsche’s  method.  As  we  will  see 
in  the  next  section,  our  convergence  estimates  scale  inversely  with  the  square  root  of 
Nitsche’s  penalty  parameter.  Hence,  we  want  to  select  Nitsche’s  penalty  parameter 
as  small  as  possible. 

6  Approximation  of  the  Generalized  Stokes  Prob¬ 
lem 

In  this  section,  we  approximate  the  homogeneous  generalized  Stokes  problem  using 
the  discrete  velocity  and  pressure  spaces  introduced  in  the  previous  section.  We  prove 
continuity,  stability,  and  a  priori  error  estimates  for  our  discretization  scheme  in  the 
single  patch  setting,  and  we  explicitly  track  all  of  our  estimates’  dependencies  on  the 
problem  parameters  and  Nitsche’s  penalty  parameter. 

6.1  Variational  Formulation 

We  begin  this  section  by  presenting  a  discrete  variational  formulation  for  the  gener¬ 
alized  Stokes  problem.  Since  members  of  Vo,/,  do  not  satisfy  homogeneous  tangential 
Dirichlet  boundary  conditions,  we  resort  to  Nitsche’s  method  [48]  to  weakly  enforce 
no-slip  boundary  conditions.  This  requires  slightly  more  regularity  on  our  problem 
data.  Specifically,  we  assume 

v  G  W1,0O(n)  :=  {w  G  :  Vw  G  L°°(H)}  . 
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This  assumption  ensures  the  trace  of  v  on  the  boundary  <912  is  well-defined.  Now,  let 
us  define  the  following  bilinear  form: 


ah (w,  v)  =  a(w,  v)  - 


[  2z/  f  ((VAv)  n)  ■  w  +  ((Vsw)  n) 
fg rh  J F  ' 


a 


pen 


hi 


w  •  V 


ds. 


(37) 

In  the  above  equation,  Cpen  >  1  denotes  a  specially  chosen  positive  penalty  constant 
which  will  be  specified  in  the  sequel.  Our  discrete  formulation  is  written  as  follows. 


(  Find  u h  G  Vo, ft  and  ph  G  Qo,h  such  that 


(G)< 


ah( Uft,  Vft)  -  b(ph,  Vft)  +  b(qh,  uh)  =  (f,  Vft)L2(n) , 

Vv/j  G  Vq .ft,  ^ft  G  Q0,ft- 


(38) 


Note  that  the  mesh-dependent  bilinear  form  given  by  (37)  has  three  additional  terms 
in  comparison  with  the  continuous  bilinear  form:  a  penalty  term,  a  consistency  term, 
and  a  stability  term.  All  three  of  these  terms  will  prove  important  in  our  proceeding 
mathematical  analysis. 

We  have  the  following  lemma  detailing  the  consistency  of  our  numerical  method 
provided  the  exact  solution  satisfies  a  reasonable  regularity  condition. 

Lemma  6.1.  Suppose  that  the  unique  weak  solution  (u ,p)  of  (IV)  satisfies  the  regu¬ 
larity  condition  u  G  H3,/2+e(12)  for  some  e  >  0.  Then: 

ah{ u,  vh)  -  b(p ,  Vft)  +  b(qh,  u)  =  (f,  vfe)L2(n)  (39) 

for  all  vh  G  V0,h  and  qh  G  Q0)h. 


Proof.  We  trivially  have 

b(qh,  u)  =  0,  Wqh  G  Q0,ft- 

Now  let  Vft  G  Vo, ft.  By  the  trace  theorem  for  fractional  Sobolev  spaces  [57],  the 
assumption  u  G  H"!'/2+<E(r2)  guarantees  that  (Vsu)  n  is  well-defined  along  <912  and 
(Vsu)  n  G  (L2(dfl))d.  Furthermore,  (Vsv^)  n  is  well-defined  along  <912  and  (Vsv/J  n  G 
(L2(<912))d.  Hence,  the  quantity  Oft(u,  vh)  is  well-defined.  Utilizing  integration  by 
parts  and  the  fact  that  u  satisfies  homogeneous  Dirichlet  boundary  conditions  and 
Vft  satisfies  homogeneous  normal  Dirichlet  boundary  conditions,  we  have 


cift(u,  vh)  -  b(p,  Vft)  =  /  (cru  -  V  ■  (2zzVAu)  +  gradp)  •  Vftdx 

Jn 


=  /  f-Vftdx 


=  (f,Vft)La 


(Q) 


where  integration  is  to  be  understood  in  the  sense  of  distributions.  This  completes 
the  proof  of  the  lemma.  □ 
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Consistency  is  the  primary  reason  that  we  employed  Nitsche’s  method  instead  of  a 
naive  penalty  method.  Nitsche’s  method  also  admits  adjoint  consistency,  and  this  will 
allow  us  to  prove  optimal  L 2  estimates  for  our  numerical  method.  This  is  in  contrast 
with  some  standard  discontinuous  Galerkin  techniques  such  as  the  Nonsymmetric 
Interior  Penalty  Galerkin  (NIPG)  method  [50].  Furthermore,  note  that  our  method 
is  consistent  for  velocity  holds  satisfying  u  G  H^'/2+e(h2)  for  arbitrary  e  >  0.  Thus, 
our  method  is  consistent  for  such  singular  problems  as  how  over  a  backward  facing 
step.  As  a  direct  result  of  consistency,  we  have  the  following  orthogonality  condition. 

Corollary  6.1.  Let  (u h,Ph)  denote  a  solution  of  ( G ),  and  suppose  that  the  unique 
weak  solution  (u ,p)  of  (W)  satisfies  the  regularity  condition  u  G  H'^2+e(G)  for  some 
e  >  0.  Then: 

ah{ u  -  nh,  vh)  -  b(p  -  ph,  vft)  +  b(qh,  u  -  uh)  =  0  (40) 

for  all  v/j  G  V0,h  and  qh  G  Q0>h. 


Our  discretization  also  enjoys  the  following  pointwise  mass  conservation  property 
which  is  a  direct  consequence  of  Proposition  5.4. 

Corollary  6.2.  Let  (u h,Ph)  denote  a  solution  of  (G).  Then: 

clivus  =  0  (41) 


We  would  like  to  note  that  in  the  event  the  viscosity  v  vanishes  for  uniformly 
positive  a,  Problem  (G)  reduces  a  compatible  discretization  of  incompressible  Darcy 
how  subject  to  a  no-penetration  boundary  condition.  This  reduction  is  contingent 
upon  the  weak  specification  of  the  no-slip  condition.  In  this  sense,  weak  boundary 
conditions  are  essential  to  the  proper  behavior  of  the  discrete  system  under  vanishing 
viscosity.  Our  proceeding  stability  and  error  analysis  extends  trivially  to  the  case 
of  vanishing  viscosity.  Furthermore,  much  like  the  solutions  of  Navier-Stokes  hows, 
the  generalized  Stokes  solution  is  characterized  by  the  presence  of  a  sharp  boundary 
layer  for  small  u.  The  weak  no-slip  condition  alleviates  the  necessity  of  highly-refined 
boundary  layer  meshes  [6,  7,  8]. 

Remark  6.1.  If  we  wish  to  impose  non-homogeneous  tangential  Dirichlet  (e.g.,  pre¬ 
scribed  slip)  boundary  conditions,  we  add  the  following  expression  to  the  right  hand 
side  of  our  discrete  formulation: 


f Bciy h)  =  Y.  f  ((V*V/0  n)  •  UBC  +  BC  ■  Vh)  ds 

Fe  rh  Jp  \  f  J 


(42) 


where  u bc  is  a  vector  function  living  on  dLl  with  prescribed  tangential  boundary  value 
and  zero  normal  boundary  value.  The  imposition  of  non-homogeneous  normal  Dirich¬ 
let  boundary  conditions  is  executed  strongly  in  the  standard  sense. 
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6.2  Continuity  and  Stability 

We  now  establish  continuity  and  stability  estimates  for  our  discrete  formulation.  Con¬ 
tinuity  and  stability,  in  conjunction  with  consistency  and  regularity,  will  guarantee 
convergence.  To  obtain  the  estimates  in  this  subsection,  we  will  need  to  call  upon 
standard  tools  from  the  discontinuous  Galerkin  community  such  as  trace  estimates, 
and  we  will  also  rely  on  the  results  presented  in  Proposition  5.2. 

Before  proceeding,  let  us  assume  throughout  the  remainder  of  the  section  that 
the  reaction  rate  o  and  the  kinematic  viscosity  v  are  constant  over  £2.  This  will 
greatly  simplify  the  presentation  of  our  mathematical  analysis.  Nonetheless,  our 
results  extend  to  the  more  general  setting  of  variable  reaction  and  viscosity.  Let  us 
define  the  following  weighted  mesh-dependent  norm: 

IMIv(ft)  :=  <ivllH(div;n)  +2Hvl2H1(n) 

+  Ml  (Vsv)n||(i2(F))d  +  2i/  nMMI \v(F)y 

Fe r,  Fevh 

(43) 

Note  that  this  is  a  proper  norm  over  our  discrete  velocity  space  due  to  the  Poincare 
inequality 

||v||Hi(n)  <  Cpoin |v|Hi(n),  Vv  e  HM)  n  H0(div;  £2)  (44) 

where  Cpoin  is  a  positive  constant  which  depends  only  on  £2.  In  fact,  by  Proposition 
5.6  and  the  Poincare  inequality,  there  exists  a  positive  constant  Cinv  independent  of 
h,  cr,  v,  and  Cpen  such  that 

IKIIv(fc)  <  °inv  (  ^||Vft||H(div;n)  +  2z/IvIh1(U)  +2uY1  ifMMI )  (45) 

V  Fe rh  llF  J 

for  all  v/,  e  Vo,/,,.  The  above  inequality  dictates  that  our  proposed  norm  acts  as 
expected  on  the  discrete  subspace  Vo, h-  That  is,  it  is  analogous  to  a  weighted  H1-norm 
coupled  with  an  appropriate  penalty  term  to  handle  tangential  boundary  conditions. 
The  use  of  a  mesh-dependent  norm  is  fairly  standard  in  the  discontinuous  Galerkin 
community.  It  is  also  standard  in  the  stabilized  methods  community.  The  use  of  a 
weighted  norm  is  motivated  by  our  desire  to  extract  error  estimates  with  an  explicit 
dependence  on  the  problem  parameters  a  and  v  as  well  as  the  penalty  constant  Cpen. 
Let  us  define  the  following  weighted  L2-norm  for  the  pressure  space 

Iklle  :=  ^4—119111.(0),  V,  6  L&i J).  (46) 

Note  that  when  u  —  0,  our  norms  reduce  to  cr-weighted  H(div)-  and  L2-norms.  Hence, 
we  recover  the  proper  norms  for  Darcy  flow  in  the  limit  of  vanishing  viscosity. 

We  have  the  following  continuity  result. 
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Lemma  6.2.  The  following  continuity  statements  hold: 


ah( w,v)  <  C'cont||w||V(fe)||v||V(/l),  Vw.  v  G  V0, ft©  nH3/2+6(0)J  (47) 

b(p,v)  <  ||p|| q || v|| V(/i)>  VpG  Lg(fi),v  G  V0, ft  ©  (h3(^)  nH:i/2+((0)) 

(48) 

where  e  >  0  is  an  arbitrary  positive  number  and  Ccont  is  a  positive  constant  indepen¬ 
dent  of  h,  a,  v ,  Cpen ,  and  e. 


Proof.  To  establish  the  first  estimate,  we  write 


Oft(w,  v)  =  a(w,  v)  —  ^  /  2v  f  (( Vsv )  n)  ■  w  +  ((Vsw)  n)  ■  v 

Ferh  'F  ' 


a 


pen 


h, 


w  •  v 


for  some  w,  v  G  V0,ft  ©  ^Hq(Q)  D  H!,/2+,E(f2)J .  We  now  bound  ah(-,  •)  term  by  term. 
To  begin,  note  immediately  that 


o(w,  v) 


L 

Ferh  Jr 


c, 


pen 


f  hp 


w  •  V 


<  l|w||v(h)||v||v(h). 


Next,  we  write 

5]  f  M(Vv)  n)  ■  w  <  2„  £  (||  wll(L2(F))dl|  (VW)  n||(L2(F))d) 


Fe  ry 


Fe  rh 


<  m  (V'v)ni 


(L2(F))P 


hp\ 


w 


(L2(F))d 


Fe  rh 


Ferh 


<  l|w||v(/olM|v(h)- 


Similarly,  we  have 

£  /  2f  ((V6w)  n)  ■  v  <  ||w||  V(/i)||v||v(/i)- 
Fe rh 

Collecting  our  bounds,  we  have 

Oft(w,v)  <  C'cont||w||v(ft)||v||V(ft) 

with  Ccont  =  3.  To  establish  the  second  continuity  result  of  the  lemma,  we  first  write 

b(p,v)  <  ||p||L2(n)||divv||L2(n). 
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The  result  is  then  a  consequence  of 

IHl2(Q)  =  (u  +  2 u)1/2  HpIIq 

and 

||divv||La(n)  <  (a  +  2z/)~1/2  ||v||v(/i). 

□ 


Now  we  seek  a  coercivity  estimate  for  the  bilinear  form  Note  that  we 

cannot  obtain  a  general  estimate  with  respect  to  the  V(/i)-norm  because  the  bilinear 
form  imposes  no  gradient  control  along  dTl.  On  the  other  hand,  (45)  suggests  that  a 
coercivity  estimate  can  be  achieved  if  we  restrict  ourselves  to  the  discrete  space  Voy. 
To  develop  estimates  which  are  independent  of  the  problem  parameters  a  and  u,  we 
further  restrict  ourselves  to  the  divergence-free  space 

V0) h  ■=  {vh  e  V0 >h  :  divv/j  =  0}  . 

To  proceed,  we  must  make  two  assumptions  regarding  the  size  of  Cpen.  First,  in  light 
of  Proposition  5.6,  we  choose  Cpen  large  enough  such  that 

II  )  n  II  f  T,2(QK))d 

cpen  >  4hKCp0inCKorn -  (L  (dK))  ,  ViF  e  JCh,  Vh  e  V0,h  (49) 

II  V*-Nh1(*T) 

where  Cpmn  is  the  Poincare  constant  associated  with  (44)  and  Cxom  is  the  positive 
constant  associated  with  the  following  Korn’s  inequality  [11]: 

MhTO)  ^  CKorn  ( ||  V"w  ||  ^2(n))dxd  +  |' ~1^d~1')  ||  W  ||(L2(an))d)  ,  Vw  6  H^fi). 

Second,  we  assume  that 

Cpen>Ah0\dQ\-1/^  (50) 

where  ho  is  the  mesh  size  of  the  coarsest  mesh  /Co  and  <90  denotes  the  area  of  dQ. 
This  second  assumption  is  necessary  as  rotation  modes  carry  zero  energy  when  o  =  0. 
Hence,  weak  boundary  conditions  are  needed  to  control  these  modes  in  rotationally 
symmetric  (or  near  rotationally  symmetric)  configurations.  As  such  configurations 
are  of  significant  engineering  interest,  we  believe  that  any  analysis  results  should 
cover  these  situations.  Note  that  a  constant  Cpen  satisfying  the  above  assumption 
need  not  depend  on  h,  a,  or  v.  Rather,  it  only  needs  to  depend  on  the  size  of  the 
domain,  the  polynomial  degree  of  the  discretization,  the  parametric  shape  regularity, 
and  global,  mesh-invariant  measures  of  the  parametric  mapping. 

We  have  the  following  lemma  governing  the  coercivity  of  our  problem. 

Lemma  6.3.  Suppose  (49)  and  (50)  are  satisfied.  Then  we  have 

ah{wh,Wh)  >  Ccoerc\\wh\\l(h),  Vwfe  G  V0.h  (51) 

where  Ccoerc  is  a  positive  constant  independent  of  h,  a,  v,  and  Cpen. 
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Proof.  Let  w/t  e  V0,ft  be  arbitrary.  We  expand 

Wfc)  =  o( wh.wh) 


E  [  2u  f((Vswft)n) 

Ferh  ' 


wh  +  ((Vswft)n)  •  wA 


a 


pen 


hi 


W ft  • 


c, 


=  ^llw^llH(div;n)  +2^||Vswh||2L(f7))dx,i  +  2z/  E  -^11 


Pen  n  |i 2 

\\(L2(F))d 


4/v  y,  /((Vs 'wft)n)-wh 
Ferh  ^ 


^ery 


(52) 


where  we  have  used  the  divergence-free  condition  on  to  obtain  ||w/j||H(div;n)  = 
||  w/t  ||  l2(o)  •  We  now  use  Cauchy- Schwarz  to  write 


2"  V 


Fe  rh 


2  hr 


a 


pen 


4/v  y,  /  ((Vswft  )  n)  ■  wfe  < 
Ferh  ^ 

(V*wk)  n|ty(F)|J  +  ^||w»||^(W)  ■ 


(53) 


Due  to  Assumption  (49)  and  the  Poincare  inequality,  we  have 


2  h 


E  7^~ II  (VSwft)nll<W))d  <  X) 


Ferh  Cpen 


KeKh 


2C 


Korn 


IwAh1(q) 


(54) 


where  Cxom  is  the  positive  constant  (only  dependent  on  the  domain  Q)  associated 
with  the  Korn’s  inequality 


(55) 


IwIh>(i!)  ^  cK„  (||v*w||^(n))JxJ  +  |an|  1,{d  ‘’llwii^^d ,  Vw  e  h‘(si) 

Inserting  (54)  into  (53)  gives 

E  [  4z/((VswA)n)  •  wh  <  2z/  E 

J F  Ferh 


F& Th' 


2  C 


1  -|wh|Li^  + 


Korn 


IHPO)  t  2/ 


{L2(F)Y  )  ’ 


and  by  inserting  the  above  inequality  into  (52)  and  employing  (55),  we  obtain 


ah(wh,wh)  >  o'||wh||iI(div;n)  + 


C 


F/  Hi 


Korn 


H1(f2) 


,  ( Cpen  _ 2 _ \  1 1 2 

+  Kr,  U  v”^7  _  IdDIVC-b  J  II w^!  (L2(F))d" 
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Invoking  Assumption  (50),  we  have 


ah{wh,wh)  >  oHwhlli^n)  + 


+  E 

Fe  rh  * 


I  (L2(F))d 


hp  <  h  <  h0. 

The  lemma  then  follows  with  Ccoerc  =  Cf^v  min  {0.25,  0.5  (C^orn)-1}  where  Cinv  is 
the  constant  associated  with  (45).  □ 


We  need  one  more  stability  estimate.  We  need  to  satisfy  the  Babuska-Brezzi  inf- 
sup  condition.  Recall  that  we  already  proved  an  inf-sup  condition  for  our  discrete 
spaces  in  Section  5.  However,  for  that  inf-sup  condition,  we  utilized  the  H1-norm  for 
the  velocity  space.  Here,  we  must  employ  the  stronger  V(/i)-norm.  To  arrive  at  an 
inf-sup  condition  for  this  stronger  norm,  we  will  proceed  by  employing  three  pow¬ 
erful  tools:  (1)  commuting  projectors,  (2)  trace  inequalities,  and  (3)  approximation 
estimates. 


Lemma  6.4.  There  exists  a  positive  constant  (5  independent  of  h,  a,  and  v  such  that 
the  following  holds:  for  every  qh  G  Qo,h,  there  exists  a  vh  G  Voy  such  that: 

divvft  =  qh  (56) 


and 


Hence, 


llwllvpi)  < 


.  (diWfc,  qh) 

ml  sup  tt - n - r — p— 

<ih£Qo,h,qh¥=o  vheV0}h  || ^h\\v(h)  Mh\\ Q 


(57) 

(58) 


Furthermore,  the  inf-sup  constant  f3  asymptotically  scales  inversely  with  the  square 
root  of  Cpen . 


Proof.  Let  qh  G  Qo,h  be  arbitrary.  Then  we  know  there  exists  a  function  v  G  Hq(£2) 
such  that  divv  =  qh  and 

2HIvHhi(q)  ^  2uP~2hh\\2L2{n) 

where  f3  is  a  positive  constant  independent  of  v.  Let  =  n^,  v.  Then,  by  Proposi¬ 
tion  5.2,  divv/j  =  divn^,  v  =  divv  =  qh  and 

^KIIhRo)  <  2iyClMn\n)  <  ZvClP^hhWl^n) 
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(59) 


where  Cu  >  0  is  a  positive  constant  independent  of  h,  a,  77,  and  Cpen.  Similarly,  we 
have 

<Tllvfc||H(div;fi)  —  2\\qh\\2L2m.  (60) 

As  v  satisfies  homogeneous  Dirichlct  boundary  conditions,  we  can  apply  the  contin¬ 
uous  trace  inequality  (see  Theorem  3.2  of  [26])  to  obtain  the  expression 

hF  Hv/i||(i,2(F))d  =  hF  ||vh  -  v||(I2(F))d 
F&vh  Ferh 

<  °tr  {h~K2bh  -  v||2l2(k)  +  |V/l  -  v|2h1(k)) 

KdKh 

where  Ctr  is  a  positive  constant  only  dependent  on  the  shape  regularity  of  the  mesh 
family  {Mh}h<ho  and  global,  mesh- invariant  measures  of  the  parametric  mapping. 
Proposition  5.5  gives 

(h^bh  ~  v||l2W  +  |vft  -  v|^i(K))  <  Clound ||v||^1(n) 

I<&K.h 

where  Cbound  is  a  positive  constant  independent  of  h,  <7,  u ,  and  Cpen.  Thus,  we  have 

2^  ~T  llv^ll(L2(jp))d  <  2  vCl^njCfrCpenfi  2  ||  Qh  ||  ^2(q)  .  (61) 

Fe  rh  F 

Combining  (45),  (59),  (60),  and  (61),  we  have 

IKIlvpp  —  Cinvp  2  (C2  +  C2oundC2rCpen)  (a  +  2v)  \\qh\\2L2(a)- 
Hence,  (57)  holds  with 


P 


r-V 2 

^  inv 


{cl  +  c2boundcl  cpen) 


-1/2 


p. 


□ 

The  inverse  dependence  of  the  inf-sup  constant  (3  on  the  square  root  of  the  penalty 
constant  Cpen  suggests  that  Cpen  should  be  chosen  as  small  as  possible  to  satisfy 
coercivity.  Indeed,  we  have  numerically  computed  the  inf-sup  constant  /3  for  a  wide 

~  _ i  /o 

range  of  values  of  Cpen  and  found  that  the  relationship  f3  <  Cpen  is  sharp  (for 
reference,  see  the  results  listed  in  Table  1).  We  would  like  to  note  that  this  makes 
intuitive  sense  as  we  lose  inf-sup  stability  entirely  if  we  enforce  the  no-slip  condition 
strongly. 

By  Lemmata  6.2,  6.3,  and  6.4  and  Brezzi’s  Theorem  [12],  we  immediately  have 
the  following  theorem  establishing  well-posedness. 
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Table  1:  Dependence  of  the  inf-sup  constant  (5  on  Nitsche’s  penalty  constant  Cpen  for 
k'  —  1,  h  —  1/16,  and  D  =  (0,  l)2. 


r 

w pen 

1 

2 

4 

8 

16 

32 

64 

P 

order 

4.24e-l 

3.89e-l 

-0.13 

3.29e-l 

-0.24 

2.50e-l 

-0.39 

1.82e-l 

-0.46 

1.30e-l 

-0.48 

9.24e-2 

-0.49 

Theorem  6.1.  Suppose  that  the  assumptions  of  Lemma  6.3  hold  true.  Then,  Problem 
(i G )  has  a  unique  weak  solution  ( uh,ph )  G  Vo,/,  x  Qo,h-  Furthermore, 


u»||vw  <  f—  P  +  2, 'C~l)-1/2  ||f||lS(n),  (62) 

IMc  <  i  (<T  +  2 vC-*n)~1/2  +  \  ||f||La(n).  (63) 

p  \  coerc  J 


We  would  like  to  note  that  all  of  the  continuity  and  stability  estimates  here  hold 
when  v  is  identically  zero  and  o  is  positive.  Hence,  we  have  a  unified  stability  analysis 
of  Stokes  flow  and  incompressible  Darcy  flow. 

Remark  6.2.  Note  that  for  the  setting  of  constant  viscosity,  one  has 

V  ■  (2zA7su)  =  z/Au.  (64) 

This  inspires  a  different  variational  formulation  than  that  presented  here  which  is 
often  the  basis  for  numerical  discretization  (see,  for  example,  [19]).  However,  these 
discretizations  (and  their  accompanying  mathematical  analysis)  are  not  extendable  to 
the  more  difficult  and  physically  relevant  setting  of  variable  viscosity,  and  they  also 
cannot  easily  accommodate  traction  boundary  conditions. 


6.3  A  Priori  Error  Estimates 

We  are  now  ready  to  derive  a  priori  error  estimates  for  our  discrete  formulation.  We 
begin  with  the  following  lemma. 

Lemma  6.5.  Let  (u ,p)  and  (u h,Ph)  denote  the  unique  solutions  of  Problems  (' W ) 
and  (G)  respectively .  Furthermore,  assume  that  u  e  iL3//2+€(fl)  for  some  e  >  0  and 
that  the  assumptions  of  Lemma  6.3  hold  true.  Then 

||u- uh||v(ft)  <  fl  +  inf  ||u-vh||v(/i)  (65) 

V  tv  coerc  /  vh£Vo,h 

and 

Wp-PhWq  <  f1  +  i)  inf  \\p-qh\\Q  +  ^J^||u-uh||vw  (66) 

V  p  J  <lh&Qo.h  fj 
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where  Ccont  is  the  continuity  constants  given  by  Lemma  6.2,  Ccoerc  is  the  coercivity 
constant  given  by  Lemma  6.3,  and  (3  is  the  inf-sup  constant  given  by  Lemma  6.f. 


Proof.  We  first  prove  (65).  We  have  that,  for  any  v/j  e  Vo th  such  that  clivv^  =  0, 


||vh  -  uh\\y(h)  <  — - a(vh  -  uh,  vh  -  uh) 


coerc 

l 


cn 


a(vh  -  u,  wh  -  uh) 


y  coerc 

<  -  u||v(h)||vft  -  Ufc||v(ft) 

coerc 

where  we  employed  the  orthogonality  given  by  Corollary  6.1  and  the  condition 

cliv  (lift  -  V/l)  =  0 
in  the  second  line  of  the  (67).  Hence,  we  can  write 

||u  -  Uft||v(fe)  <  inf  (||u  -  Vh\\v(h)  +  ||vft  -  U/Jlypo) 

v>,£V 0,h 

<( 1  +  inf  ||u-vA||v(h). 

V  C' coerc J  vheVo ,h 

We  now  prove  (66).  We  have  that,  for  any  qt  €  Qo,h, 

1  b(ph  —  qh,  wh) 

\\Ph-qh\\Q<^  sup 


(67) 


(68) 


(3  w,,eVo,ft  ||wA||V(h) 

1  b(p  -  qh,  wft)  -  a(u  -  uh,  wh) 

=  —  SUP  - 77 - 77 - 

(3  whev0 ,h  l|Wh||v(h) 

<  i  {\\P~  Qh\\Q  +  Ccont\\u-uh\\v{h)) 

h' 


(69) 


where  we  again  employed  orthogonality  in  the  second  line  above.  Inequality  (66)  then 
follows  in  the  same  manner  as  (68)  by  a  splitting  of  the  pressure  error  and  a  usage  of 
(69).  □ 


We  have  the  following  theorem  giving  us  a  priori  convergence  estimates  which  are 
optimal  for  the  discrete  velocity  field  and  suboptimal,  by  one  order,  for  the  discrete 
pressure  field. 

Theorem  6.2.  Let  (u ,p)  and  (u h,Ph)  denote  the  unique  solutions  of  Problems  (W) 
and  ( G )  respectively.  Furthermore,  assume  that  (u ,p)  e  HJ+1(f2)  x  7P(f2)  for  some 
j  >  1/2  and  that  the  assumptions  of  Lemma  6.3  hold  true.  Then 

||u  -  uft||v(h)  <  cu  (l  +  ^F!L\  Vcrh2s+ 2  +  2z//r2s||u||HS+i(n)  (70) 

\  ^ coerc J 


31 


and 


\\p  -  Ph\\<2  <  Cp  (l  +  Vj  (a  +  2v)  1/2  hs\\p\\Hs(n)  +  ^^\\u-uh\\v{h)  (71) 

for  s  =  min  where  k'  is  the  polynomial  degree  of  our  discretization,  Cu  is  a 

positive  constant  independent  of  h,  a,  and  v  which  asymptotically  scales  with  the 
square  root  of  Cpen,  and  Cp  is  a  positive  constant  independent  of  h,  a,  v,  and  Cpen. 


Proof.  We  first  prove  (70).  Recall  the  error  estimate  given  by  (65): 

Ilu  -  uh||v(h)  <  (l  +  inf  ||u  -  vh||V(h). 

V  C^coercJ  vheV0,d 


Noting  divll^u  =  IIgh  clivu  =  0,  we  can  choose  =  Ily  u  in  the  above  expression 
to  obtain 


u 


Ufc||v(fc)  <  1  + 


C, 


cont 


cn 


u 


nvhullv(ft) 


=  Cco^T{  +  T2  +  T3  +  r4  (72) 

where  we  have  assigned  Cco  =  ( 1  +  ffont  )  and 

V  h^coerc  J 


T\  = 

=  °\\' 

u  — 

nVhUllH(div;Q) 

=  u||u-  n°hu||^2(n) 

(73) 

t2  = 

=  2u 

u  — 

TT°  nl2 

11vauIh1(q) 

(74) 

t3  = 

=  2u 

E 

hF ||  (Vs  (u- 

nvh  U))  nll(L2(F))d 

(75) 

F&h 

T,= 

=  2v 

E 

Cpen  hp  1 1 U 

nVhUll(I,2(F))d. 

(76) 

Ferh 


To  handle  the  face  integral  in  (75),  we  recruit  the  multiplicative  trace  inequality 
for  fractional  Sobolev  spaces  [57]  and  Young’s  inequality  element-wise  to  obtain  the 
bound 


Y  CpenhF ||  (Vs  (u  -  n°ft  u))  n\\2{L2{F))d  < 

Ferh 

(Ctrc,l)  (jU  —  ny^tlln1^)  +  h^\u  —  nyfeU|H9+i^^ 

K&Kh 

where  1/2  <  q  <  s  and  Ctrc,i  is  a  positive  constant  independent  of  h,  a,  u,  and 
Cpen.  To  handle  the  face  integral  in  (76),  we  recruit  the  standard  continuous  trace 
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inequality  element-wise  to  obtain  the  bound 


'y  ^  Cpenhp  ||  U  IIyh  u||  < 

F&Vh 

(Ctrc, 2?  (,lJ<  HU  “  nV,UllL2(K)  +  lU  -  nV,UlHhK)) 

KeKh 

where  Ctrc,2  is  a  positive  constant  independent  of  h,  cr,  and  v  which  varies  linearly 
with  the  square  root  of  Cpen.  It  should  be  noted  the  two  constants  Ctrc,  1  and  Ctrc,2 
necessarily  depend  on  the  shape  regularity  of  the  mesh  family  {Q}h<ho  and  the  para¬ 
metric  mapping  which  together  give  the  shape  regularity  of  the  mesh  family  {/C}^<^0 . 
See  [26]  for  more  details.  Inserting  the  above  two  inequalities  into  (72)  and  then 
applying  Proposition  5.5,  we  immediately  acquire  the  bound 

||u  -  n£hu||v(h)  <  CuCcoV ah2s+2  +  2i//i2s||u||hh-i(q) 

for  Cu  a  positive  constant  independent  of  h,  cr,  and  v  with  the  same  functional 
dependency  on  the  penalty  parameter  as  Ctrc, 2- 

The  proof  for  (71)  is  much  more  immediate.  Choosing  qh  =  II Qhp  in  the  error 
estimate  given  by  (66),  one  obtains 

\\P-Ph\\<2  <  (i  +  i)  lb  -  n^p||e  +  ^j±\\u-uh\\V(h)- 

Inequality  (71)  follows  by  an  application  of  Proposition  5.5  to  bound  the  pressure 
interpolation  error.  □ 

Since  we  have  the  bound 


2^|  •  IhPq)  II  '  ||v(/i)j 

the  above  theorem  also  provides  optimal  convergence  rates  for  the  velocity  field  in 
the  H1  -norm.  If  we  assume  slightly  more  regularity  for  the  pressure  space,  we  have 
the  following  proposition. 

Proposition  6.1.  Let  (u ,p)  and  (u h,Ph)  denote  the  unique  solutions  of  Problems 
(W)  and  (G)  respectively.  Furthermore,  assume  that  (u ,p)  e  HJ+1  (O)  x  7P+1  (Q) 
for  some  j  >  1/2  and  that  the  assumptions  of  Lemma  6.3  hold  true.  Then 

\\p-Ph\\Q  <  Cp,e  ^1  +  {a  +  2vYl/2  hs+l\\p\\Ha+i{n)  +  -  uh\\v{h)  (77) 

for  s  =  min  {k' ,  j}  where  Cp,e  is  a  positive  constant  independent  of  h,  a,  v,  and  Cpen. 
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Observe  that  the  pressure  error  estimate  given  by  the  above  proposition  is  still 
suboptimal  clue  to  the  presence  of  the  velocity  error,  which  converges  with  order  s  for 
general  viscous  flows.  Let  us  further  note  that  the  preceding  theorem  and  proposition 
are  trivially  extended  to  the  setting  of  vanishing  viscosity.  In  this  case,  the  velocity 
error  actually  converges  with  order  s  +  1,  giving  optimal  a  priori  error  estimates  for 
both  the  discrete  pressure  held  and  discrete  velocity  held  for  incompressible  Darcy 
how. 

Under  an  elliptic  regularity  assumption,  we  can  obtain  optimal  estimates  for  the 
velocity  held  in  the  L2-norm  by  utilizing  a  standard  duality  argument.  Given  the 
unique  solutions  (u ,p)  and  (uh,ph)  of  Problems  (W)  and  (G),  let  us  consider  the 
following  ancillary  problem,  written  in  strong  form. 

Find  (xp,r)  £  Hq(D)  x  L'l(Q)  such  that 

axp  —  V  •  (2v\7sxp)  +  gradr  =  u  —  uh  in  Q 

divxp  =  0  in  Q 

xp  =  0  on  d hi 

The  above  problem  has  a  unique  weak  solution  (xp,r).  Before  proceeding,  note  that 
we  can  formally  take  the  divergence  of  (78)  to  obtain 

A r  =  0,  in  Q.  (81) 

Since  r  has  zero  average,  it  follows,  at  least  from  our  formal  argument,  that  r  =  0. 
This  argument  can  be  made  rigorous  by  a  suitable  use  of  convolutions  and  passing 
to  the  limit.  Now  suppose  that  xp  G  H2(D).  We  can  then  multiply  the  left  and  right 
hand  sides  of  (78)  by  a xp  +  V  •  (2vWsxp)  to  acquire  the  result 

WIIl2^)  +  II V  •  (2^VSV>)  1 1  l2  (q)  =  (u  -  uh,  axp  +  V  ■  (2z/V^))l2(0)  •  (82) 

A  simple  application  of  Cauchy-Schwarz  and  the  triangle  inequality  gives 

H^llL2(n)  +  llv  ■  (2^VA'i/>)  ||22(n)  = 
l|u  -  ufc||L2(n)  (lk^llL2(n)  +  || V  ■  (2z/VAh/0  ||L2(n))  (83) 

Since  x2+y2  >  \(x  +  y)2)  we  can  divide  both  sides  by  + 1| \7 •  (2zv Vs-*/?) 

to  obtain 

\  ( WlliAn)  +  llv  '  (2zA7AV>)  || l2(st2) )  <  llu  -  w||l2(Q)-  (84) 

As  xp  satishes  normal  and  tangential  homogeneous  Dirichlet  boundary  conditions 
and  v  is  assumed  positive,  we  can  employ  a  combination  of  Korn’s  inequalities  and 
Poincare  inequalities  to  obtain  a  standard  elliptic  regularity  result  of  the  form 

ll^llH2(a)  <  0_1||u  -  uh||La(n).  (85) 

where  Ca  is  a  positive  constant  which  only  depends  on  the  domain  Q.  In  view  of  the 
above  discussion,  we  have  the  following  theorem. 


(78) 

(79) 

(80) 
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Theorem  6.3.  Let  (u ,p)  and  (u h,Ph)  denote  the  unique  solutions  of  Problems  ( W ) 
and  ( G )  respectively,  and  let  (i f,r)  denote  the  unique  solution  of  Problem  ( A ).  Fur¬ 
thermore,  assume  that  (u,p)  G  H?+1  (Q)  x  7P(f2)  for  some  j  >  1,  that  if  G  H2(h2)? 
and  that  the  assumptions  of  Lemma  6.3  hold  true.  Then 

llu  -  w||l2(Q)  <  C^s+1||u||  HS+1(T2)  (86) 

for  s  =  min  {k',  j}  where  Ci  is  a  positive  constant  independent  of  h,  a,  and  v  which 
asymptotically  scales  with  the  square  root  of  Cpen. 


Proof.  Since  by  assumption  ( if,r )  G  HJ(f2)  x  i/1(f2),  consistency  and  symmetry  give 

ah(v,if)  —  b(r ,  v)  =  (u  -  uh,  v)La(n) 
for  all  v  G  Vft.  Let  us  take  v  =  u  —  u^.  We  then  have 

II  u  -  uj22(n)  =  ah(  u  -  uh,  if). 

as  div  (u  —  Uh)  =  0  (or  r  =  0  by  our  preceding  discussion).  By  using  the  orthogonality 
given  by  Corollary  6.1,  we  can  write 

llu  -  u4l2(Q)  =  a^(u  -  Uhi'lf  -  Kh^) 

—  Ccont  ||u-U/l||V(h)||V?-  nvh^||v(h)-  (87) 

We  bound  the  interpolation  error  by  utilizing  a  similar  argument  to  that  used  to 
prove  (77),  obtaining 

H  -  uvh^\\  v(h)  <  Cinterp  ( a1/2h 2  +  vl/2h)  1 1 V7 1 1  h2  («) 

for  Cinterp  a  positive  constant  independent  of  h,  a,  and  v  which  asymptotically  scales 
with  the  the  square  root  of  Cpen.  We  can  now  employ  the  elliptic  regularity  condition 
(85)  to  arrive  at 

|| if  ~  n^^||v(h)  <  CACinterp  (ct1 1 2 V~X if1  +  U~1/2h)  ||u  -  U/,||L2(n).  (88) 

Inserting  (88)  into  (87)  results  in 

||  U  U/j, 1 1  l2(Q)  ^  C AC inter pC cont  (&  ^  V  T  O  ^  tl)  ||u  U/j||y(/q. 

Immediately  invoking  the  error  estimate  given  by  Theorem  6.2,  we  obtain 

II u  U/, || 5;  Ctemp  fl  G  y/Dah)  h  1 1 u 1 1  jjs+1  (q)  (89) 
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where  CtemV  is  a  positive  constant  independent  of  h,  a,  and  v  which  asymptotically 
scales  with  the  square  root  of  Cpen  and 


Dah  =  — 
v 

However,  we  also  have  the  following  estimate  due  to  Theorem  6.2: 

IIu-u^IIl2^)  <C*u  (1  +  ^^)  (1  +  (JDa/l)_1)1/2hs+1||u||HS+i(n).  (90) 

\  ^ coerc  / 

The  desired  result  follows  by  taking  the  minimum  of  (89)  and  (90).  □ 


This  concludes  our  a  priori  error  analysis.  Note  that,  for  reasonably  regular 
exact  solutions,  we  have  obtained  optimal  estimates  for  the  velocity  held  in  both 
the  strong  V(/i)-norm  as  well  as  the  weaker  H1-  and  L2-norms.  The  estimates  are 
additionally  robust  with  respect  to  the  fluid  coefficients  a  and  v.  On  the  other 
hand,  we  have  obtained  pressure  error  estimates  which  are  suboptimal  by  one  order. 
This  is  reminiscent  of  error  estimates  for  stabilized  equal-order  interpolations  of  the 
Stokes  equations  and  is  not  unexpected  as  both  our  discrete  velocity  and  pressure 
spaces  consist  of  mapped  piecewise  polynomials  which  are  only  complete  up  to  degree 
k' .  However,  our  later  numerical  studies  suggest  the  conservative  nature  of  these 
estimates  by  revealing,  for  simple  model  problems,  optimal  convergence  rates  for  the 
pressure  held.  Ongoing  work  is  being  dedicated  to  the  theoretical  confirmation  of 
these  convergence  rates.  Note  that  our  analysis  covers  typical  singular  solutions  of 
the  generalized  Stokes  equations.  Later  in  this  chapter,  we  will  numerically  study  the 
effectiveness  of  our  method  for  a  selection  of  singular  Stokes  problems.  Finally,  we 
would  like  to  mention  that  our  velocity  error  estimates  are  completely  independent 
of  the  pressure  held.  This  property  does  not  hold  for  discretizations  which  preserve 
the  incompressibility  in  only  a  discrete  sense. 

Remark  6.3.  In  opposition  to  standard  Bubnov- Galerkin  methods,  the  constant 

(l  +  Ccont  ( Ccoerc )  ) 

appearing  in  (65)  cannot  be  reduced  to  just  Ccont  (C'coerc)_1. 

Remark  6.4.  At  this  point,  suitable  requirements  guaranteeing  elliptic  regularity  for 
domains  obtained  by  NURBS  mappings  are  unknown.  One  anticipates  that  such  re¬ 
quirements  should  be  less  stringent  than  those  associated  with  polyhedral  domains 
(namely,  convexity)  because  of  enhanced  smoothness.  In  our  view,  this  is  an  interest¬ 
ing  area  of  research. 
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Figure  8:  Example  multi-patch  construction  in  M2. 

7  Extension  to  Multi-Patch  Domains 

As  was  mentioned  previously  in  Subsection  4.3,  most  geometries  of  scientific  and 
engineering  interest  cannot  be  represented  by  a  single  patch.  Instead,  the  multi¬ 
patch  concept  must  be  invoked.  We  assume  that  there  exist  np  sufficiently  smooth 
parametric  mappings  Fj  :  (0,  l)d  — y  such  that  the  subdomains 

Oj  Fj  t  1)  •  •  •  j  Tip 

are  non-overlapping  and 

n  =  uyjTi. 

We  refer  to  each  subdomain  Q,:  (and  its  pre-image)  as  a  patch.  For  a  visual  depiction 
of  a  multi-patch  construction  in  M2,  see  Figure  8.  We  build  discrete  velocity  and 
pressure  spaces  over  each  patch  Qj,  i  =  1  ,...,np  in  the  same  manner  as  in  the 
previous  sections  except  that  we  do  not  yet  enforce  boundary  conditions,  and  we 
denote  these  spaces  as  W(f2j)  and  Qh(fb). 

To  proceed  further,  we  must  make  some  assumptions.  First  of  all,  we  assume  that 
if  two  disjoint  patches  f2j  and  Qj  have  the  property  that  d fij  D  dflj  ^  0,  then  this 
intersection  consists  strictly  of  patch  faces,  edges,  and  corners.  More  succinctly,  two 
patches  cannot  intersect  along  an  isolated  portion  of  a  face  (or  edge)  interior.  Second, 
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we  assume  that  the  mappings  {F,;}^1  are  compatible  in  the  following  sense:  if  two 
patches  ftj  and  ftj  share  a  face,  then  F,  and  F;  parametrize  that  face  identically  up 
to  changes  in  orientation.  Third,  we  assume  that  if  two  patches  ftj  and  fl,-  share  a 
face,  the  B-spline  meshes  associated  with  the  patches  are  identical  along  that  face. 
This  guarantees  our  mesh  is  conforming.  Finally,  we  assume  for  simplicity  that  k\  = 
. . .  =  kd  =  k*  for  all  patches.  The  mixed  polynomial  degree  case  introduces  additional 
complications  that  are  beyond  the  scope  of  this  work.  We  would  like  to  note  that  all 
four  assumptions  hold  if  we  employ  a  conforming  NURBS  multi-patch  construction. 
See,  for  example,  Chapter  2  of  [20]. 

We  define  our  global  discrete  velocity  and  pressure  spaces  as  follows: 

V0,h  ■■=  {vh  G  H0  (div;  ft)  :  vja.  G  14(ft?:),  Vi  =  1, ... ,  np}  ,  (91) 

Qo,h  '■=  {<lh  £  Lo  (ft)  :  Qh\ni  £  Qh(^i),  Vi  =  1, . . .  ,npj  .  (92) 

The  space  Vo y  is  easily  constructed  due  to  our  preceding  four  assumptions  and  use 
of  open  knot  vectors.  Specifically,  we  set  to  zero  the  coefficient  of  any  basis  function 
whose  normal  is  nonzero  along  <9ft,  and  along  shared  faces  between  patches,  we  (i) 
equivalence  the  coefficients  of  any  basis  functions  whose  normal  values  are  nonzero 
and  equal  in  magnitude  and  direction  and  (ii)  set  opposite  the  coefficients  of  any 
basis  functions  whose  normal  values  are  nonzero,  equal  in  magnitude,  and  opposite 
in  direction.  We  note  that  this  is  precisely  the  same  procedure  as  is  used  to  construct 
Raviart-Thomas  spaces  on  conforming  finite  element  meshes.  We  simply  have  patches 
instead  of  elements.  It  is  easily  shown  that  the  spaces  Vo y  and  along  with  the 
divergence  operator,  form  the  bounded  discrete  cochain  complex 

)  div 

Uo  ,h  - >  Wlo,/!- 

However,  functions  in  Vo y  do  not  necessarily  he  in  H1(ft)  as  tangential  continuity  is 
not  enforced  across  patch  interfaces.  Hence,  we  need  to  account  for  this  lack  of  conti¬ 
nuity  when  designing  a  discretization  scheme  for  the  generalized  Stokes  equations.  We 
employ  the  symmetric  interior  penalty  method  [1,  24,  62],  a  standard  technique  in  the 
discontinuous  Galerkin  community,  to  weakly  enforce  tangential  continuity  between 
adjacent  patches. 

We  now  establish  some  preliminary  notation.  Let  /C/j  (ft*)  and  V-)<(ftj)  denote  the 
sets  of  physical  mesh  elements  and  faces  associated  with  patch  ftj.  We  denote  the 
global  set  of  mesh  elements  as  fCh  and  the  global  set  of  mesh  faces  as  Fh-  As  in  the 
single  patch  setting,  we  define  the  boundary  mesh  to  be 

Th  =  {F  G  ^(ftj),?  =  1, .  , .  ,np  :  F  C  <9ft}  ,  (93) 

and  we  define  the  interface  mesh  to  be 

lh  =  {F  G  JX(fti),  i  =  1  ,...,np:Fe  ±  j  and  F  £Yh}  .  (94) 
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For  each  face  F  G  X/-(  belonging  to  the  interface  mesh,  there  exist  two  unique  adjacent 
elements  K+,  K~  G  ICh  such  that  F  G  dK+  and  F  G  dK~ .  We  define  for  such  a  face 
the  mesh  size 

Iif  ~  ( hK+  +  hK~) .  (95) 

Let  (j)  be  an  arbitrary  scalar-,  vector-,  or  matrix-valued  piecewise  smooth  function, 
and  let  us  denote  by  (j)+  and  (j)~  the  traces  of  (j)  on  F  as  taken  from  within  the  interior 
of  K+  and  K~  respectively.  We  define  the  mean  value  of  0  at  x  G  F  as 

M  :=  \  (0+  +  <t>~ )  •  (96) 

Further,  for  a  generic  multiplication  operator  0,  we  define  the  jump  of  0©n  at  x  G  F 
as 

[</>  0  n]  :=  (j)+  0  nA'+  +  (\T  0  nK-  (97) 

where  11^-+/-  denotes  the  outward  facing  normal  on  the  boundary  d K+/~  of  element 
K+/~. 

With  the  above  notation  established,  let  us  define  the  following  bilinear  form: 


i=  1 


<(w>v)  =  J]  ({2uVsw,  Vsv)(L2(Q.)}dxd  +  (^w,v)l2(q^ 

L  ' 

[  2u  (§Vsv}}  :  [w  0  n]  +  §Vsw}  :  [v  0  n])  ds 


Fexh 


Feih  JF 

£  /  2x((fVv)n) 
fg rh  Jf  k 


-FF  [w  ®  n]  :  [v  0  nl  ]  ds 

hp 


a 


■  w  +  ((Vsw)  n)  •  v  —  "  w  •  v  )  ds.  (98) 


Above,  Cpen  >  0  denotes  the  same  positive  penalty  constant  as  before.  Our  discrete 
formulation  over  the  multi-patch  domain  then  reads  as  follows. 

(  Find  11^  G  Voy,  and  ph  G  such  that 


(MP)  { 


a*h{uh,vh)  -b(ph,vh)  +  b{qh,  ufc)  =  (f ,v/l)L2(n)  (99) 


[  for  all  \h  G  V0l h  and  qh  G  Q0)h. 


As  in  the  single  patch  setting,  the  discrete  formulation  detailed  above  returns  a  point- 
wise  divergence-free  velocity  field.  However,  we  do  not  have  a  convergence  analysis 
available  as  we  do  not  yet  have  a  multi-patch  analogue  of  Proposition  5.2.  We  antic¬ 
ipate  this  will  take  new  theoretical  developments.  Nonetheless,  we  have  utilized  the 
above  formulation  in  practice  and  observed  it  returns  optimal  convergence  rates  for 
both  velocity  and  pressure  fields. 
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8  Numerical  Verification  of  Convergence  Estimates 

In  this  section,  we  numerically  verify  our  convergence  estimates  using  a  collection  of 
problems  with  exact  solutions.  Throughout,  we  choose  Nitsche’s  penalty  constant  as 

Cpen  —  5  (k‘  +  1) 

where  k'  is  the  polynomial  degree  of  a  given  discretization.  We  have  found  that  this 
choice  leads  to  stable  numerical  formulations  for  the  generalized  Stokes  equations. 
Furthermore,  unless  otherwise  specified,  we  employ  uniform  parametric  meshes,  linear 
parametric  mappings,  and  B-spline  spaces  of  maximal  continuity. 

8.1  Two-Dimensional  Manufactured  Solution 

As  a  first  numerical  experiment,  we  consider  a  two-dimensional  manufactured  solution 
that  was  originally  presented  in  [14],  Let 

Q  =  (0,1)2 


and 

f  =  au  -  V  ■  (2zA7su)  +  Vp 

with 

__  2ex(— 1  +  x)2x2(y2  —  y)(— 1  +  2y) 

U  (— ex(— 1  +  x)x(—2  +  x(3  +  x))(— 1  +  y)2y2) 

and 

p  =  (-424  +  156e  +  (c/2-y)(-456  +  ea;(456  +  a;2(228  -  5(c/2-c/))  + 

2cc(— 228  +  ( y 2  -  y))  +  2a;3(-36  +  ( y 2  -  y))  +  a;4(12  +  ( y 2  -  y))))). 

Homogeneous  boundary  conditions  are  applied  along  the  boundary  <912,  and  the  pres¬ 
sure  is  enforced  to  satisfy  fnpdx  =  0.  The  unique  solution  to  the  generalized  Stokes 
equation  with  the  prescribed  forcing  is  then  clearly  (u,p)  =  (u ,p).  The  streamlines 
and  pressure  contours  associated  with  the  exact  solution  are  plotted  in  Figure  9.  Note 
from  the  streamline  plot  that  the  velocity  solution  has  a  simple  vortex  structure. 

For  the  above  constructed  solution,  we  have  computed  convergence  rates  for 
divergence-conforming  B-splinc  discretizations  of  varying  mesh  size  and  polynomial 
degree.  Furthermore,  we  have  computed  convergence  rates  for  a  variety  of  Damkohler 
numbers 


v 


where  L  is  a  length  parameter  which  we  henceforth  specify  as  one.  These  convergence 
rates  are  provided  in  Tables  2,  3,  and  4.  Note  immediately  from  the  tables  that  our 
theoretically  derived  error  estimates  are  confirmed.  Second,  note  that  the  L2-norm 
of  the  pressure  error  optimally  converges  like  0(hk+1),  which  is  an  improvement 
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X 


(a) 


(b) 

Figure  9:  Generalized  Stokes  manufactured  solution  in  2-D:  (a)  Flow  velocity  stream¬ 
lines,  (b)  Pressure  contours. 
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Table  2:  Generalized  Stokes  convergence  rates  in  2-D:  Da  =  0 


Polynomial  degree  k'  —  1 


h 

1/4 

1/8 

1/16 

1/32 

1/64 

u  —  uh||v(h) 

7.75e-2 

3.96e-2 

1.98e-2 

9.90e-3 

4.95e-3 

order 

- 

0.97 

1.00 

1.00 

1.00 

u  —  uh  HpO) 

5.48e-2 

2.80e-2 

1.40e-2 

7.00e-3 

3.50e-3 

order 

- 

0.97 

1.00 

1.00 

1.00 

u  —  u/i  L2(Q) 

2.77e-3 

8.16e-4 

2.28e-4 

6.10e-5 

1.58e-5 

order 

- 

1.76 

1.84 

1.90 

1.95 

\\P'-Ph\\iP(fl) 

5.04e-3 

1.38e-3 

3.49e-4 

8.72e-5 

2.18e-5 

order 

- 

1.87 

1.98 

2.00 

2.00 

Polynomial  degree  k'  —  2 


h 

1/4 

1/8 

1/16 

1/32 

1/64 

||u  -  Uh||v(h) 

1.37e-2 

3.30e-3 

8.03e-4 

1.98e-4 

4.92e-5 

order 

- 

2.06 

2.04 

2.02 

2.01 

u  —  Uh  H^n) 

9.70e-3 

2.33e-3 

5.68e-4 

1.40e-4 

3.48e-5 

order 

- 

2.06 

2.04 

2.02 

2.01 

u  —  uh  L2(Q) 

2.94e-4 

3.84e-5 

5.03e-6 

6.47e-7 

8.21e-8 

order 

- 

2.94 

2.93 

2.96 

2.98 

\\P  -  Ph\\v*(fl) 

1.18e-3 

1.19e-4 

1.17e-5 

1.19e-6 

1.27e-7 

order 

- 

3.31 

3.35 

3.30 

3.23 

Polynomial  degree  k'  —  3 


h 

1/4 

1/8 

1/16 

1/32 

1/64 

||u  —  Uft||v(h) 

1.39e-3 

1.81e-4 

2.33e-5 

3. Ole-6 

3.85e-7 

order 

- 

2.94 

2.96 

2.95 

2.97 

u  —  uh  HpO) 

9.83e-4 

1.28e-4 

1.65e-5 

2.10e-6 

2.66e-7 

order 

- 

2.94 

2.96 

2.97 

2.98 

u  —  L2(Q) 

3.05e-5 

2.34e-6 

1.59e-7 

1.03e-8 

6.55e-10 

order 

- 

3.70 

3.88 

3.95 

3.98 

\\P  ~  Ph\\L*{0) 

1.10e-4 

5.64e-6 

3.45e-7 

2.19e-8 

1.39e-9 

order 

- 

4.29 

4.03 

3.98 

3.98 
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Table  3:  Generalized  Stokes  convergence  rates  in  2-D:  Da  =  1 


Polynomial  degree  k'  —  1 


h 

1/4 

1/8 

1/16 

1/32 

1/64 

u  —  uh||v(h) 

7.75e-2 

3.96e-2 

1.98e-2 

9.90e-3 

4.95e-3 

order 

- 

0.97 

1.00 

1.00 

1.00 

u  —  uh  HpO) 

5.47e-2 

2.80e-2 

1.40e-2 

7.00e-3 

3.50e-3 

order 

- 

0.97 

1.00 

1.00 

1.00 

u  —  u/i  L2(Q) 

2.75e-3 

8.09e-4 

2.26e-4 

6.05e-5 

1.57e-5 

order 

- 

1.76 

1.84 

1.90 

1.95 

\\P'-Ph\\iP(fl) 

5.04e-3 

1.37e-3 

3.48e-4 

8.72e-5 

2.18e-5 

order 

- 

1.88 

1.98 

2.00 

2.00 

Polynomial  degree  k'  —  2 


h 

1/4 

1/8 

1/16 

1/32 

1/64 

||u  -  Uh||v(h) 

1.37e-2 

3.30e-3 

8.03e-4 

1.98e-4 

4.92e-5 

order 

- 

2.06 

2.04 

2.02 

2.01 

u  —  Uh  H^n) 

9.70e-3 

2.33e-3 

5.68e-4 

1.40e-4 

3.48e-5 

order 

- 

2.06 

2.04 

2.02 

2.01 

u  —  uh  L2(Q) 

2.94e-4 

3.84e-5 

5.03e-6 

6.47e-7 

8.21e-8 

order 

- 

2.94 

2.93 

2.96 

2.98 

\\P  -  Ph\\v*(fl) 

1.18e-3 

1.19e-4 

1.17e-5 

1.19e-6 

1.27e-7 

order 

- 

3.31 

3.35 

3.30 

3.23 

Polynomial  degree  k'  —  3 


h 

1/4 

1/8 

1/16 

1/32 

1/64 

||u  —  Uft||v(h) 

1.39e-3 

1.81e-4 

2.35e-5 

3. Ole-6 

3.85e-7 

order 

- 

2.94 

2.95 

2.96 

2.97 

u  —  uh  HpO) 

9.83e-4 

1.28e-4 

1.65e-5 

2.10e-6 

2.66e-7 

order 

- 

2.94 

2.96 

2.97 

2.98 

u  —  L2(Q) 

3.05e-5 

2.34e-6 

1.59e-7 

1.03e-8 

6.55e-10 

order 

- 

3.70 

3.88 

3.95 

3.98 

\\P  ~  Ph\\L*{0) 

1.10e-4 

5.64e-6 

3.45e-7 

2.19e-8 

1.39e-9 

order 

- 

4.29 

4.03 

3.98 

3.98 
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Table  4:  Generalized  Stokes  convergence  rates  in  2-D:  Da  =  1000 


Polynomial  degree  k'  —  1 


h 

1/4 

1/8 

1/16 

1/32 

1/64 

u  —  uh||v(h) 

3.13e-3 

1.33e-3 

6.36e-4 

3.14e-4 

1.57e-4 

order 

- 

1.23 

1.06 

1.02 

1.00 

u  —  uh  HpO) 

5.50e-2 

2.78e-2 

1.39e-2 

6.98e-3 

3.49-3 

order 

- 

0.98 

1.00 

0.99 

1.00 

u  —  u/i  L2(Q) 

1.93e-3 

4.79e-4 

1.22e-4 

3.20e-5 

8.38e-6 

order 

- 

2.01 

1.97 

1.93 

1.93 

\\P'-Ph\\iP(fl) 

3.37e-3 

7.97e-4 

1.96e-4 

4.89e-5 

1.22e-5 

order 

- 

2.08 

2.02 

2.00 

2.00 

Polynomial  degree  k'  —  2 


h 

1/4 

1/8 

1/16 

1/32 

1/64 

||u  -  Uh||v(h) 

5.24e-4 

l.lle-4 

2.58e-5 

6.30e-6 

1.55e-6 

order 

- 

2.24 

2.11 

2.03 

2.02 

u  —  Uh  H^n) 

9.85e-3 

2.32e-3 

5.67e-4 

1.40e-4 

3.48e-5 

order 

- 

2.08 

2.03 

2.02 

2.01 

u  —  uh  L2(Q) 

2.83e-4 

3.80e-5 

5. Ole-6 

6.46e-7 

8.21e-8 

order 

- 

2.90 

2.92 

2.96 

2.98 

\\P  -  Ph\\v*(fi) 

4.69e-4 

5.36e-5 

6.42e-6 

7.97e-7 

9.98e-8 

order 

- 

3.13 

3.06 

3.01 

3.00 

Polynomial  degree  k'  —  3 


h 

1/4 

1/8 

1/16 

1/32 

1/64 

||u  —  Uft||v(h) 

5.31e-5 

6.12e-6 

7.57e-7 

9.56e-8 

1.22e-8 

order 

- 

3.12 

3.02 

2.99 

2.97 

u  —  uh  HpO) 

9.74e-4 

1.26e-4 

1.64e-5 

2.10e-6 

2.66e-7 

order 

- 

2.95 

2.94 

2.97 

2.98 

u  —  L2(Q) 

3.04e-5 

2.33e-6 

1.59e-7 

1.03e-8 

6.55e-10 

order 

- 

3.71 

3.89 

3.95 

3.98 

\\p~Ph\\v(n) 

7.56e-5 

4.60e-6 

3.19e-7 

2.12e-8 

1.37e-9 

order 

- 

4.04 

3.89 

3.91 

3.95 
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Figure  10:  Generalized  Stokes  manufactured  solution  in  3-D:  Flow  velocity  stream¬ 
lines  colored  by  velocity  magnitude. 

over  our  theoretically  derived  estimate.  Third,  note  that  our  methodology  is  robust 
with  respect  to  the  Damkohler  number.  That  is,  the  errors  for  our  discretization 
are  virtually  independent  of  the  Damkohler  number.  In  fact,  our  pressure  error 
decreases  with  increasing  Damkohler  number.  Finally,  it  should  be  mentioned  that 
(a)  the  H1  error  of  the  velocity  field  approaches  the  H1  best  approximation  error 
as  k'  is  increased,  and  (b)  the  L2  error  of  the  pressure  held  approaches  the  L2  best 
approximation  error  as  k'  is  increased. 

8.2  Three-Dimensional  Manufactured  Solution 

As  a  second  numerical  experiment,  we  consider  a  three-dimensional  manufactured 
solution  representing  a  vortical  filament.  Let 

n  =  (o,i)3 


and 

with 


f  =  cru  -  V  ■  (2iA7"u)  +  Vp 


0  = 


u  =  curl</>, 

x(x  —  1  )y2(y  —  1  )2z2(z  —  l)2 

0 

x2(x  —  l)2y2(y  —  1  )2z(z  —  1) 
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Table  5:  Generalized  Stokes  convergence  rates  in  3-D:  Da  =  0 


Polynomial  degree  k'  —  1 


h 

1/2 

1/4 

1/8 

1/16 

1/32 

u  —  uh\ \v(h) 

2.59e-2 

1.27e-2 

5.91e-3 

2.81e-3 

1.36e-3 

order 

- 

1.03 

1.10 

1.07 

1.05 

U  —  U/j)  H1(f2) 

1.51e-2 

7.64e-3 

3.77e-3 

1.87e-3 

9.34e-4 

order 

- 

0.98 

1.02 

1.01 

1.00 

u  —  uh  L2(Q) 

1.35e-3 

3.68-4 

1.03e-4 

2.81e-5 

7.40e-6 

order 

- 

1.88 

1.84 

1.87 

1.93 

\\v~Vh\W{n) 

5.41e-2 

1.48e-2 

3.58e-3 

8.85e-4 

2.26e-4 

order 

- 

1.87 

2.05 

2.02 

1.97 

Polynomial  degree  k'  —  2 


h 

1/2 

1/4 

1/8 

1/16 

1/32 

||u  -  Uh||v(h) 

6.50e-3 

1.54e-3 

4.10e-4 

9.51e-5 

2.15e-5 

order 

- 

2.08 

1.91 

2.11 

2.15 

u  —  uh  HpO) 

3.71e-3 

9.90e-4 

2.79e-4 

6.59e-5 

1.50e-5 

order 

- 

1.91 

1.83 

2.08 

2.14 

u  —  u7i  II  L2(Q) 

1.97e-4 

4.25e-5 

7.38e-6 

8.67e-7 

9.18e-8 

order 

- 

2.21 

2.53 

3.09 

3.23 

\\p-Phh^Q) 

1.50e-2 

1.59e-3 

2.00e-4 

2.56e-5 

3.26e-6 

order 

- 

3.24 

2.99 

2.97 

2.97 

and 

4 

p  =  sin(7rx)  sin(7n/) - -. 

7 H 

Again,  homogeneous  boundary  conditions  are  applied  along  the  boundary  dfl,  and 
the  pressure  is  enforced  to  satisfy  fnpdx  =  0.  The  unique  solution  to  the  generalized 
Stokes  equation  with  the  prescribed  forcing  is  then  (u,p)  =  (u ,p).  Streamlines  asso¬ 
ciated  with  the  exact  solution  are  plotted  in  Figure  10.  Note  that  the  the  streamlines 
wrap  around  a  single  diagonal  vortex  filament. 

As  in  the  two-dimensional  setting,  we  have  computed  convergence  rates  for  a 
variety  of  divergence-conforming  B-spline  discretizations  and  Damkohler  numbers 


Da 


crL2 

v 


where  L  is  a  length  parameter  which  we  again  specify  as  equal  to  one.  These  conver¬ 
gence  rates  are  summarized  in  Tables  5,  6,  and  7.  Note  immediately  from  the  tables 
that  our  theoretically  derived  error  estimates  are  confirmed.  Second,  note  that  the 
L2-norm  of  the  pressure  error  optimally  converges  like  0(hk  +1),  which  is  an  improve¬ 
ment  over  our  theoretically  derived  estimate.  Third,  note  that  our  method  is  robust 
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Table  6:  Generalized  Stokes  convergence  rates  in  3-D:  Da  =  1 


Polynomial  degree  k'  —  1 


h 

1/2 

1/4 

1/8 

1/16 

1/32 

||u  -  Ufc||v(/») 

2.59e-2 

1.28e-2 

5.91e-3 

2.81e-3 

1.36e-3 

order 

- 

1.02 

1.11 

1.07 

1.05 

u  —  Uh  Hpn) 

1.51e-2 

7.64e-3 

3.77e-3 

1.87e-3 

9.33e-4 

order 

- 

0.98 

1.02 

1.01 

1.00 

u  —  U/1  L2(Q) 

1.34e-3 

3.66e-4 

1.02e-4 

2.79e-5 

7.34e-6 

order 

- 

1.87 

1.84 

1.87 

1.93 

\\P  ~  Ph\\iP(n) 

5.41e-2 

1.48e-2 

3.58e-3 

8.85e-4 

2.21e-4 

order 

- 

1.87 

2.05 

2.02 

2.00 

Polynomial  degree  k'  —  2 


h 

1/2 

1/4 

1/8 

1/16 

1/32 

u  —  Uh||v(h) 

6.50e-3 

1.54e-3 

4.10e-4 

9.50e-5 

2.15e-5 

order 

- 

2.08 

1.91 

2.11 

2.14 

U  —  uh  HpO) 

3.71e-3 

9.89e-4 

2.79e-4 

6.59e-5 

1.50e-5 

order 

- 

1.91 

1.83 

2.08 

2.14 

u  —  L2(Q) 

1.97e-4 

4.24e-5 

7.38e-6 

8.64e-7 

9.18e-8 

order 

- 

2.22 

2.52 

3.09 

3.23 

\\p-Ph\\L*(n) 

1.50e-2 

1.59e-3 

2.00e-4 

2.56e-5 

3.26e-6 

order 

- 

3.24 

2.99 

2.97 

2.97 
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Table  7:  Generalized  Stokes  convergence  rates  in  3-D:  Da  =  1000 


Polynomial  degree  k'  —  1 


h 

1/2 

1/4 

1/8 

1/16 

1/32 

||u  -  Uh||v(h) 

1.55e-3 

5.16e-4 

2.04e-4 

9.07e-5 

4.32e-5 

order 

- 

1.59 

1.34 

1.17 

1.07 

u  —  Uh  Hpn) 

1.56e-2 

7.57e-3 

3.74e-3 

1.86e-3 

9.31e-4 

order 

- 

1.04 

1.02 

1.01 

1.00 

u  —  U/1  L2(Q) 

1.16e-3 

2.58e-4 

6.29e-5 

1.59e-5 

4.11e-6 

order 

- 

2.17 

2.04 

1.98 

1.95 

\\P  ~  Ph\\iP(n) 

5.41e-2 

1.48e-2 

3.57e-3 

8.84e-4 

2.20e-4 

order 

- 

1.87 

2.05 

2.01 

2.01 

Polynomial  degree  k'  —  2 


h 

1/2 

1/4 

1/8 

1/16 

1/32 

u  —  Uh||v(h) 

2.98e-4 

5.25e-5 

1.06e-5 

2.42e-6 

5.84e-7 

order 

- 

2.50 

2.31 

2.13 

2.05 

U  —  uh  HpO) 

3.79e-3 

8.56e-4 

2.08e-4 

5.14e-5 

1.28e-5 

order 

- 

2.15 

2.04 

2.02 

2.01 

u  —  L2(Q) 

1.88e-4 

2.84-5 

3.74e-6 

4.97e-7 

6.06e-8 

order 

- 

2.73 

2.93 

2.91 

3.04 

\\p-Ph\\L*(n) 

1.50e-2 

1.59e-3 

2.00e-4 

2.56e-5 

3.26e-6 

order 

- 

3.24 

2.99 

2.97 

2.97 
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with  respect  to  the  Damkohler  number.  Finally,  let  us  remark  the  exact  velocity  field 
is  recovered  for  discretizations  of  degree  k'  >  3. 

8.3  Two-Dimensional  Problem  with  a  Singular  Solution 

To  examine  how  our  discretization  performs  in  the  presence  of  singularities,  we  con¬ 
sider  Stokes  flow  in  the  L-shaped  domain  Yt  =  ( — 1, 1)2)\([0, 1)  x  (—1,0]).  The  flow 
problem  in  consideration  is  depicted  in  Figure  11(a).  Homogeneous  Dirichlet  bound¬ 
ary  conditions  are  applied  along  T d  =  { (0,  y)  :  y  G  (—1,0)}  U  {(#,  0)  :  x  G  (0,1)}, 
Neumann  boundary  conditions  are  applied  along  Y m  —  <9Q\r p,  and  we  set  a  =  0, 
v  —  1,  and  f  =  0.  As  in  [60],  the  Neumann  boundary  conditions  are  chosen  such  that 
the  exact  solution  is 

_  rA((l  +  A)  sin (6)^(6)  +  cos (9)^' (9) 

11  rA(— (1  +  A)  cos (d)'ip(d)  +  sin 

and 

p  =  — rA_1  ((1  +  X)2(p'(6)  +  (f>"'(9))  /(I  —  A) 

where  ( r,6 )  are  polar  coordinates  with  respect  to  the  origin  (0,0), 

cj)  =  sin((l  +  \)6)  cos(Au;)/(l  +  A)  —  cos((l  +  \)6) 

—  sin((l  —  A )9)  cos(Acu)/(l  —  A)  +  cos((l  —  A )d), 

u>  =  !</>,  and  A  ~  0.54448373678246  is  the  smallest  positive  root  of 

sin(Ao;)  +  Asin(o;)  =  0. 

This  singular  solution  is  illustrated  in  Figure  12.  Note  that  (u,p)  G  H1+A(0)  x 
Hx(Yl).  This  is  the  strongest  corner  singularity  for  the  Stokes  operator  in  the  L- 
shaped  domain,  and,  as  such,  this  numerical  example  models  typical  singular  behavior 
observed  in  the  vicinity  of  reentrant  corners. 

To  compute  this  flow  example  using  our  discretization  technique,  we  must  resort 
to  a  multi-patch  construction.  We  utilize  the  three-patch  construction  illustrated  in 
Figure  11(b).  Each  patch  is  mapped  from  the  parametric  domain  using  an  affine 
parametrization.  As  discussed  in  Section  7,  we  impose  normal  continuity  strongly  be¬ 
tween  patches  and  tangential  continuity  weakly  using  the  symmetric  interior  penalty 
method.  We  have  computed  convergence  rates  for  a  variety  of  divergence-conforming 
B-spline  discretizations  and  reported  our  results  in  Table  8.  Note  from  the  table  that 
the  V(/i)-norm  of  the  velocity  field  and  the  L2-norm  of  the  pressure  field  are  approach¬ 
ing  the  optimal  convergence  rates  of  0(hx )  as  h  — >  0.  Furthermore,  note  the  velocity 
and  pressure  errors  improve  with  increasing  polynomial  degree.  This  is  somewhat 
counterintuitive  as  we  also  increase  smoothness  with  polynomial  degree.  This  being 
said,  such  a  property  has  also  been  observed  in  the  context  of  Maxwell’s  equations 
[15].  While  we  only  employed  uniform  meshes  for  the  computations  reported  here, 
one  could  of  course  obtain  more  satisfactory  results  with  geometrically  graded  meshes 
[30,  31], 
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(b) 

Figure  11:  Singular  Stokes  solution  in  2-D:  (a)  Problem  setup,  (b)  Multi-patch  con¬ 
struction. 
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(b) 

Figure  12:  Singular  Stokes  solution  in  2-D:  (a)  Flow  velocity  streamlines,  (b)  Pressure 
contours. 
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Table  8:  Singular  Stokes  convergence  rates  in  2-D 


Polynomial  degree  k'  —  1 


h 

1/2 

1/4 

1/8 

1/16 

1/32 

u  —  u/i| \v(h) 

2.29e0 

1.61e0 

1.13e0 

7.73e-l 

5.33e-l 

order 

- 

0.51 

0.51 

0.55 

0.54 

\\p  -  Ph\\v*(n) 

1.30e0 

9.40e-l 

6.62e-l 

4.66e-l 

3.32e-l 

order 

- 

0.47 

0.51 

0.51 

0.49 

Polynomial  degree  k'  —  2 


h 

1/2 

1/4 

1/8 

1/16 

1/32 

u  —  u/i|  V(/i) 

1.50e0 

1.05e0 

7.25e-l 

5.00e-l 

3.46e-l 

order 

- 

0.51 

0.53 

0.54 

0.53 

\\P  -  PhWvvi) 

8.56e-l 

6.39e-l 

4.50e-l 

3.20e-l 

2.36e-l 

order 

- 

0.42 

0.51 

0.49 

0.44 

Polynomial  degree  k'  —  3 


h 

1/2 

1/4 

1/8 

1/16 

1/32 

|  U  —  Ufo||v(/i) 

1.15e0 

8.36e-l 

5.79e-l 

4.00e-l 

2.78e-l 

order 

- 

0.46 

0.53 

0.53 

0.52 

\\p-Ph\\LHO) 

6.05e-l 

4.90e-l 

3.54e-l 

2.57e-l 

1.91e-l 

order 

- 

0.30 

0.47 

0.46 

0.43 
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8.4  Cylindrical  Couette  Flow 


Couette  flow  is  often  used  as  a  “sanity  check”  for  Stokes  and  Navier-Stokes  discretiza¬ 
tions.  Cylindrical  Couette  flow  is  a  more  realistic  problem  which  describes  the  flow 
between  two  concentric  rotating  cylinders.  Here,  we  consider  generalized  Stokes  flow 
between  a  fixed  outer  cylinder  and  a  rotating  inner  cylinder.  The  problem  setup  is 
illustrated  in  Figure  13.  No  external  forcing  is  applied.  In  the  absence  of  Darcy  drag 
forces  (be.,  a  =  0),  the  velocity  field  for  this  flow  assumes  the  form 

ue(r)  sm(0) 

U  ue(r)  cos(9) 


where 

,  v  ,  B 

ue(r)  —  Ar -\ - , 

r 

(r,  9)  are  polar  coordinates  with  respect  to  the  center  of  the  cylinders,  and 


A  =  -a 


<52 


1  -<52’ 


b  =  n. 


(1  —  <52)  ’ 


O  -IL  A  -  r> 

1  L{n  j  0 


T'out 


We  have  depicted  this  velocity  field  in  Figure  14(a).  In  the  presence  of  Darcy  drag 
forces,  the  character  of  the  flow  field  changes  considerably.  Notably,  the  motion  of 
the  fluid  is  confined  to  a  small  boundary  layer  attached  to  the  inner  cylinder.  This 
motion  explicitly  takes  the  form 


we(r)sin(0) 

Ug  (v )  COs(d) 

where 

=  jj  h  (Tr )K\  (jr out)  -  h  (jrout)Ki  ( jr ) 
h  (7 rin)K1  (yrout)  -  h  (jrout)K1  (7 rin)  ’ 

7  =  \J o ~Jv,  and  I\  and  K\  are  modified  Bessel  functions  of  the  first  and  second  kind 
respectively.  Note  that  7_1  acts  as  a  length  scale,  and  the  width  of  the  boundary  layer 
attached  to  the  inner  cylinder  is  proportional  to  7-1.  We  have  depicted  a  velocity 
field  corresponding  to  7  =  y/50  in  Figure  14(b).  Finally,  under  the  constraint  that 


0, 


the  pressure  field  is  identically  zero  for  both  the  Stokes  limit  as  well  as  the  generalized 
Stokes  setting.  In  what  follows,  we  assume  rm  =  1,  rout  =  2,  and  U  —  1. 

We  have  computed  convergence  rates  for  a  variety  of  divergence-conforming  B- 
spline  discretizations  and  for  7  =  0  and  7  =  y/bO.  To  represent  the  annular  domain 
in  our  computations,  we  employed  the  polar  mapping 


F(6,6) 


{(rout  -  Dn)6  +  rin)  sin^TT^) 
((rout  -  Dn)6  +  rin)  cos(27t£i) 


,V(6,6)e(o,i)2 


(100) 
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out 

Figure  13:  Cylindrical  Couette  flow:  Problem  setup. 


and  periodic  B-splines  of  maximal  continuity  in  the  ^-direction  (see  Section  2  of  [25]). 
It  should  be  emphasized  that  we  do  not  use  the  polar  form  of  the  generalized  Stokes 
equations.  Rather,  we  utilize  the  polar  mapping  to  define  our  divergence-conforming 
B-splines  in  physical  space  and  then  employ  the  Cartesian-based  variational  formu¬ 
lation  discussed  in  this  chapter.  The  results  of  our  computations  are  summarized 
in  Tables  9  and  10.  Note  from  the  tables  that  all  of  onr  theoretically  derived  error 
estimates  are  confirmed,  though  the  results  corresponding  to  7  =  a/50  have  a  more 
substantial  pre-asymptotic  range  due  to  the  presence  of  a  boundary  layer.  Addition¬ 
ally,  note  that  we  obtain  null  pressure  fields  and  axisymmetric  velocity  fields  with 
null  radial  component. 

We  repeated  our  computations  using  the  multi-patch  NURBS  construction  illus¬ 
trated  in  Figure  16.  Each  of  the  four  patches  are  built  through  a  sufficient  rotation 
of  the  canonical  quadratic  single-element  NURBS  patch  described  in  Figure  17,  and 
we  have  tabulated  the  location  and  weights  of  the  control  points  of  the  canonical 
quadratic  patch  in  Figure  11.  The  resulting  NURBS  parametrization  is  identical  to 
the  polar  parametrization  in  the  radial  direction  and  different  in  the  angular  direction. 
We  define  our  B-spline  discretization  scheme  on  the  multi-patch  NURBS  construc¬ 
tion  using  the  procedure  outlined  in  Section  7  with  normal  continuity  being  enforced 
strongly  between  patches  and  tangential  continuity  being  enforced  weakly.  Surpris¬ 
ingly,  we  found  we  obtained  exactly  the  same  velocity  and  pressure  fields  using  the 
multi-patch  NURBS  construction  as  we  did  with  the  polar  mapping.  This  property  is 
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Figure  14:  Cylindrical  Couette  flow:  (a)  Flow  velocity  arrows  for  7  =  0,  (b)  Flow 
velocity  arrows  for  7  =  \/50- 
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h/ho  =  1/8 

Figure  15:  Cylindrical  Couette  flow:  Sequence  of  polar  meshes. 
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Table  9:  Cylindrical  Couette  flow  convergence  rates:  7  =  0 


Polynomial  degree  k'  —  1 


h/h0 

1/2 

1/4 

1/8 

1/16 

1/32 

||U  -  Uh||v(h) 

5.42e-l 

2.63e-l 

1.26e-l 

6.12e-2 

2.99e-2 

order 

- 

1.04 

1.06 

1.04 

1.03 

u  —  uh  HflO) 

4.48e-l 

2.32e-l 

1.17e-l 

5.86e-2 

2.93e-2 

order 

- 

0.95 

0.99 

1.00 

1.00 

u  —  u/i  L2(Q) 

5.00e-2 

1.53-2 

4.28e-3 

1.14e-3 

2.94e-4 

order 

- 

1.71 

1.84 

1.91 

1.96 

\\ur  —  (ur)h\\L2(a) 

0 

0 

0 

0 

0 

\\P  -  Ph\\vnn) 

0 

0 

0 

0 

0 

Polynomial  degree  k'  —  2 


h/h0 

1/2 

1/4 

1/8 

1/16 

1/32 

||u  -  Uh||v(h) 

9.77e-2 

2.42e-2 

5.64e-3 

1.32e-3 

3.14e-4 

order 

- 

2.01 

2.10 

2.10 

2.07 

u  —  Uh  Hfln) 

7.68e-2 

2.00e-2 

4.92e-3 

1.21e-3 

2.99e-4 

order 

- 

1.94 

2.02 

2.02 

2.02 

u  —  uh  L2(Q) 

4.43e-3 

6.03e-4 

8.13e-5 

1.07e-5 

1.38e-6 

order 

- 

2.88 

2.89 

2.93 

2.95 

\\ur  ~  (ur)h\ \l2(Q.) 

0 

0 

0 

0 

0 

\\p  -  Ph\\v*m 

0 

0 

0 

0 

0 

Polynomial  degree  k'  —  3 


h/ho 

1/2 

1/4 

1/8 

1/16 

1/32 

||u  -  Ufc||v(/») 

2. Ole-2 

2.67e-3 

3.27e-4 

4. Ole-5 

4.98e-6 

order 

- 

2.91 

3.03 

3.03 

3.01 

u  —  uh  Hfln) 

1.52e-2 

2.13e-3 

2.84e-4 

3.72e-5 

4.80e-6 

order 

- 

2.84 

2.91 

2.93 

2.95 

u  —  uh  L2(Q) 

6.59e-4 

5.69e-5 

4.82e-6 

3.50e-7 

2.33e-8 

order 

- 

3.53 

3.56 

3.78 

3.91 

|| i^r)h  ||l/2(0) 

0 

0 

0 

0 

0 

\\p-Ph\\mn) 

0 

0 

0 

0 

0 
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Table  10:  Cylindrical  Couette  flow  convergence  rates:  7  = 


Polynomial  degree  k'  —  1 


h/ho 

1/2 

1/4 

1/8 

1/16 

1/32 

u  —  u/i| \v(h) 

8.06e-l 

4.59e-l 

2.30e-l 

1.10e-l 

5.28e-2 

order 

- 

0.81 

1.00 

1.06 

1.06 

u  —  uh  HpO) 

3.75e0 

2.44e0 

1.35e0 

6.93e-l 

3.49e-l 

order 

- 

0.62 

0.85 

0.96 

0.99 

u  —  uh  L2(Q) 

2.43e-l 

9.90e-2 

3.16e-2 

9.24e-3 

2.56e-3 

order 

- 

1.30 

1.65 

1.77 

1.85 

| \ur  ~  (ur)h\\L2(Q) 

0 

0 

0 

0 

0 

\\P  ~  Ph\\u>(Sl) 

0 

0 

0 

0 

0 

Polynomial  degree  k'  —  2 


h/ho 

1/2 

1/4 

1/8 

1/16 

1/32 

||u  -  Uh||v(h) 

3.97e-l 

1.30e-l 

3.28e-2 

7.56e-3 

1.74e-3 

order 

- 

1.61 

1.99 

2.12 

2.12 

u  ~~  uh  HpO) 

1.96e0 

6.99e-l 

1.86e-l 

4.57e-2 

1.12e-2 

order 

- 

1.49 

1.91 

2.03 

2.03 

u  —  L2(Q) 

9.86e-2 

2.00e-2 

2.73e-3 

3.66e-4 

4.86e-5 

order 

- 

2.30 

2.87 

2.90 

2.91 

|| {^r)h  ||l2(0) 

0 

0 

0 

0 

0 

\\p-Ph\\mn) 

0 

0 

0 

0 

0 

Polynomial  degree  k'  —  3 


h/ho 

1/2 

1/4 

1/8 

1/16 

1/32 

||u  -  Uh||v(h) 

1.69e-l 

3.15e-2 

4.  Ole-3 

4.80e-4 

5.86e-5 

order 

- 

2.42 

2.97 

3.06 

3.03 

u  —  Uh  Hpn) 

8.54e-l 

1.67e-l 

2.25e-2 

2.95e-3 

3.86e-4 

order 

- 

2.35 

2.89 

2.93 

2.93 

u  —  U/1  L2(Q) 

3.31e-2 

3.56e-3 

3.03e-4 

2.55e-5 

1.83e-6 

order 

- 

3.22 

3.55 

3.57 

3.80 

||^r  iPjr)h  ||l2(0) 

0 

0 

0 

0 

0 

\\p~Ph\\LHn) 

0 

0 

0 

0 

0 

Patch  2 


Figure  16:  Cylindrical  Couette  flow:  NURBS  multi-patch  construction 
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Figure  17:  Cylindrical  Couette  flow:  Patch  template  control  points 
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Tabic  11:  Cylindrical  Couette  flow:  Patch  template  control  points. 


Control  Point 

X 

y 

w 

1 

0 

1 

1 

2 

1 

1 

1/V2 

3 

1 

0 

1 

4 

0 

3/2 

1 

5 

3/2 

3/2 

1/V2 

6 

3/2 

0 

1 

7 

0 

2 

1 

8 

2 

2 

1/V2 

9 

2 

0 

1 

purely  a  consequence  of  the  Piola  transform.  To  shed  some  light  on  this  observation, 
let  us  consider  a  parametric  velocity  field  of  the  form  v  =  {u(£2),0}T.  Then,  the 
image  of  this  velocity  field  under  the  Piola  transform,  denoted  as  v,  has  the  property 
that 

1 - ^ 

divv  =  —  divv  =  0 
u 

where  J  is  the  determinant  of  the  Jacobian  matrix  D F.  Moreover,  the  vector  v  is 
oriented  in  the  direction  of  parametric  lines  defined  by  £2  =  C  where  C  is  an  arbitrary 
constant.  Hence,  if  represents  the  angular  direction  and  £2  represents  the  radial 
direction,  this  implies  that  divergence-free  lateral  velocity  fields  in  parametric  space 
are  mapped  to  axisymmetric  angular  velocity  fields  in  physical  space  regardless  of 
the  specification  of  F.  Unfortunately,  this  argument  applies  only  to  angular  velocity 
fields  and  we  are  not  generally  able  to  obtain  axisymmetric  pressure  and  radial  velocity 
fields  using  the  multi-patch  NURBS  construction. 


8.5  Annular  Poiseuille  Flow 


Poiseuille  flow  is  another  problem  which  is  often  utilized  as  a  “sanity  check”  for  Stokes 
and  Navier-Stokes  discretizations.  Here,  we  consider  generalized  Stokes  Poiseuille  flow 
of  a  viscous  fluid  between  two  concentric  cylinders.  The  problem  setup  is  illustrated 
in  Figure  18.  No-slip  and  no-penetration  boundary  conditions  are  imposed  along 
the  cylinder  surfaces,  and  periodic  boundary  conditons  are  imposed  along  the  axial 
direction.  An  external  axial  pressure  gradient  is  applied  to  drive  the  fluid.  In  the 
absence  of  Darcy  drag  forces,  the  velocity  held  for  annular  Poiseuille  how  assumes 
the  form 


u  = 


0 

0 


uz(r) 
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Figure  18:  Annular  Poiseuillc  flow:  Problem  setup. 


61 


where 


r2  _  2 

1  out  ' in 


Uz  r  = 


Ap 


2  2  i 

r  ■  —  r  4- 

'in  1 


hi(r  out/'r  *- 


ln  (r/ry 


(r,  0)  are  polar  coordinates  with  respect  to  the  center  of  the  cylinders,  and  —  Ap/L  is 
the  applied  pressure  gradient.  In  the  presence  of  Darcy  drag,  the  velocity  held  takes 
the  form 


u  = 


0 

0 


uz(r) 


where 


A p  f  Kpjjr)  ( logout )  ~  jo(7Dn))  ~  IqM  {K0(jrout)  -  Ko(jrin)) 
crL  V  K0(jrout)I0(jrin)  -  K0(jrin)I0(rrout) 

A  p 
crL’ 


7  =  \fa~Jv,  and  J0  and  K0  are  modihed  Bessel  functions  of  the  hrst  and  second  kind 
respectively.  We  have  illustrated  velocity  fields  corresponding  to  7  =  0  and  7  =  VbO 
in  Figure  19.  As  was  the  case  for  cylindrical  Couette  how,  7-1  acts  as  a  length  scale, 
and  there  is  a  boundary  layer  attached  to  the  inner  cylinder  with  width  proportional 
to  7-1.  Finally,  under  the  constraint  that 


/  pdx  =  0, 

Jn 

the  pressure  held  is  identically  zero  for  both  the  Stokes  limit  as  well  as  the  generalized 
Stokes  setting.  In  what  follows,  we  assume  rtn  =  1,  ront  =  2,  L  =  1,  and  A p  =  —1. 

We  have  computed  convergence  rates  for  a  variety  of  divergence-conforming  B- 
spline  discretizations  and  for  7  =  0  and  7  =  We  have  employed  both  the  polar 
mapping  and  the  NURBS  multi-patch  construction  described  in  the  previous  subsec¬ 
tion  to  represent  the  annular  domain  in  our  computations  with  the  axial  direction 
parametrized  using  a  simple  linear  mapping.  Periodic  B-splines  of  maximal  continu¬ 
ity  were  employed  in  the  axial  direction.  The  results  of  our  computations  for  7  =  0 
are  summarized  in  Tables  12  and  13.  Note  immediately  that  all  of  our  theoretically 
derived  error  estimates  are  confirmed.  Additionally,  note  that  we  obtain  null  pressure 
helds,  radial  velocity  helds,  and  angular  velocity  holds  for  all  cases  considered.  This 
being  said,  observe  that  we  obtain  slightly  smaller  axial  velocity  errors  with  the  polar 
mapping  than  we  do  with  the  NURBS  multi-patch  construction.  This  is  because  we 
maintain  exact  axisymmetry  with  the  polar  mapping  but  not  with  the  NURBS  multi¬ 
patch  construction.  The  results  of  our  computations  for  7  =  a/50  are  summarized 
in  Tables  14  and  15.  Again,  note  that  our  theoretically  derived  error  estimates  are 
confirmed,  and  our  results  for  the  polar  mapping  induce  slightly  smaller  axial  velocity 
errors  than  our  results  for  the  NURBS  multi-patch  construction. 
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(b) 


Figure  19:  Annular  Poiseuille  flow:  Axial  velocity  contours  along  an  axial  slice  for 
(a)  7  =  0  and  (b)  7  =  \/50. 
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Table  12:  Annular  Poiseuillc  flow  convergence  rates  for  7  =  0:  Polar  meshes 


Polynomial  degree  k'  —  1 


h/h0 

1/2 

1/4 

1/8 

1/16 

1/32 

u  —  Uh||v(h) 

5.30e-l 

2.47e-l 

1.19e-l 

5.80e-2 

2.86e-2 

order 

- 

1.10 

1.05 

1.04 

1.02 

u  —  uh  HpO) 

4.56e-l 

2.28e-l 

1.14e-l 

5.67e-2 

2.83e-2 

order 

- 

1.00 

1.00 

1.01 

1.00 

1  u  —  u/i  L2(Q) 

5.49e-2 

1.62-2 

4.34e-3 

1.12e-3 

2.84e-4 

order 

- 

1.76 

1.90 

1.95 

1.98 

|| Hr  i^r)h  ||l2(0) 

0 

0 

0 

0 

0 

ll^e  —  iue)h\ \l2(Q) 

0 

0 

0 

0 

0 

\\P  ~  PhWtfiSl) 

0 

0 

0 

0 

0 

Polynomial  degree  k'  —  2 


h/h0 

1/2 

1/4 

1/8 

1/16 

1/32 

u  —  u/i|  V(/i) 

5.67e-2 

1.23e-2 

2.83e-3 

6.73e-4 

1.64e-4 

order 

- 

2.20 

2.12 

2.07 

2.04 

U  —  uh  HpO) 

4.49e-2 

1.07e-2 

2.61e-3 

6.45e-4 

1.60e-4 

order 

- 

2.07 

2.04 

2.02 

2.01 

u  —  u/i|  L2(Q) 

2.41e-3 

3.49e-4 

4.64e-5 

5.97e-6 

7.58e-7 

order 

- 

2.79 

2.91 

2.96 

2.98 

\\ur  —  {ur)h\  Ll2(0) 

0 

0 

0 

0 

0 

\\ue  ~  (ue)h\ U2(H) 

0 

0 

0 

0 

0 

\\p-Ph\\LHQ) 

0 

0 

0 

0 

0 

Polynomial  degree  k'  —  3 


h/h0 

1/2 

1/4 

1/8 

1/16 

1/32 

||u  -  Uh||v(h) 

1.86e-3 

2.14e-4 

2.54e-5 

3.11e-6 

3.87e-7 

order 

- 

3.12 

3.07 

3.03 

3.01 

u  —  Uh  H^n) 

1.37e-3 

1.74e-4 

2.26e-5 

2.93e-6 

3.75e-7 

order 

- 

2.98 

2.94 

2.95 

2.97 

u  —  u/i|  L2(Q) 

5.38e-5 

5.05e-6 

4.03e-7 

2.80e-8 

1.83e-9 

order 

- 

3.41 

3.65 

3.85 

3.94 

\\ur  {^r)h  ||l2(Q) 

0 

0 

0 

0 

0 

\\U0  ~  (ue)h  |  l2(Q) 

0 

0 

0 

0 

0 

\\p-Ph\\mn) 

0 

0 

0 

0 

0 
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Table  13:  Annular  Poiseuillc  flow  convergence  rates  for  7  =  0:  NURBS 


Polynomial  degree  k'  —  1 


h/ho 

1/2 

1/4 

1/8 

1/16 

1/32 

u  —  Uh||v(h) 

5.30e-l 

2.47e-l 

1.19e-l 

5.80e-2 

2.86e-2 

order 

- 

1.10 

1.05 

1.04 

1.02 

u  —  H  1{1Ch) 

4.57e-l 

2.28e-l 

1.14e-l 

5.67e-2 

2.83e-2 

order 

- 

1.00 

1.00 

1.01 

1.00 

1  u  —  u/i  L2(Q) 

5.51e-2 

1.63-2 

4.35e-3 

1.12e-3 

2.84e-4 

order 

- 

1.76 

1.91 

1.96 

1.98 

|| Hr  i^r)h  ||l2(0) 

0 

0 

0 

0 

0 

ll^e  —  iue)h\ \l2(Q) 

0 

0 

0 

0 

0 

\\P  ~  PhWtfiSl) 

0 

0 

0 

0 

0 

Polynomial  degree  k'  —  2 


h/ho 

1/2 

1/4 

1/8 

1/16 

1/32 

u  —  u/i|  V(/i) 

5.69e-2 

1.23e-2 

2.83e-3 

6.75e-4 

1.64e-4 

order 

- 

2.21 

2.12 

2.07 

2.04 

u  -  uh  Hp/QO 

4.51e-2 

1.07e-2 

2.61e-3 

6.46e-4 

1.61e-4 

order 

- 

2.08 

2.04 

2.01 

2.00 

u  —  u/i|  L2(Q) 

2.52e-3 

3.55e-4 

4.70e-5 

6.04e-6 

7.66e-7 

order 

- 

2.83 

2.92 

2.96 

2.98 

\\ur  —  (' Ur)h\\L2(Q ) 

0 

0 

0 

0 

0 

\\ue  ~  (ue)h\ U2(H) 

0 

0 

0 

0 

0 

\\p-Ph\\LHQ) 

0 

0 

0 

0 

0 

Polynomial  degree  k'  —  3 


h/ho 

1/2 

1/4 

1/8 

1/16 

1/32 

||u  -  uh||v(/o 

1.91e-3 

2.38e-4 

2.73e-5 

3.33e-6 

4.13e-7 

order 

- 

3.00 

3.12 

3.04 

3.01 

u  _  H  1(JCh) 

1.44e-3 

2.03e-4 

2.48e-5 

3.16e-6 

4.03e-7 

order 

- 

2.83 

3.03 

2.97 

2.97 

u  —  uh|  L2(Q) 

8.00e-5 

1.21e-5 

6.41e-7 

4.00e-8 

2.53e-9 

order 

- 

2.72 

4.24 

4.00 

3.98 

\\ur  {^r)h  ||l2(Q) 

0 

0 

0 

0 

0 

\\ue  —  (ue)h  |  L2(Q) 

0 

0 

0 

0 

0 

\\p-Ph\\mn) 

0 

0 

0 

0 

0 
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Table  14:  Annular  Poiseuillc  flow  convergence  rates  for  7  =  Polar  meshes 


Polynomial  degree  k'  —  1 


h/h0 

1/2 

1/4 

1/8 

1/16 

1/32 

||u-uh||V(h) 

1.73e-l 

9.88e-2 

4.91e-2 

2.35e-2 

1.13e-2 

order 

- 

0.81 

1.01 

1.06 

1.06 

u  —  Ufc  h1^) 

8.20e-l 

5.29e-l 

2.90e-l 

1.48e-l 

7.46e-2 

order 

- 

0.65 

0.87 

0.97 

0.99 

u  —  u^l  L2(Q) 

4.88e-2 

2.12-2 

6.69e-3 

1.94e-3 

5.37e-4 

order 

- 

1.20 

1.66 

1.79 

1.85 

\\ur  —  (' Ur)h\\L2(Cl ) 

0 

0 

0 

0 

0 

Me  —  (ue)h\ U2(fi) 

0 

0 

0 

0 

0 

||p-Pfc|U2(n) 

0 

0 

0 

0 

0 

Polynomial  degree  k'  —  2 


h/h0 

1/2 

1/4 

1/8 

1/16 

1/32 

|  U  —  Uh||v(h) 

8.77e-2 

2.59e-2 

6.37e-3 

1.46e-3 

3.37e-4 

order 

- 

1.76 

2.02 

2.13 

2.12 

u  —  uh  HpO) 

4.34e-l 

1.39e-l 

3.63e-2 

8.86e-3 

2.16e-3 

order 

- 

1.64 

1.94 

2.03 

2.04 

u  —  uh  L2(Q) 

2.32e-2 

3.98e-3 

5.33e-4 

7. 14e-5 

9.46e-6 

order 

- 

2.54 

2.90 

2.90 

2.92 

|| (Ur)h  ||l2(0) 

0 

0 

0 

0 

0 

\We  ~  (ue)h\ \l2(Q) 

0 

0 

0 

0 

0 

\\p-Ph\\l?(Sl) 

0 

0 

0 

0 

0 

Polynomial  degree  k'  —  3 


h/ho 

1/2 

1/4 

1/8 

1/16 

1/32 

u  —  Uh||v(h) 

2.95e-2 

5.72e-3 

7.10e-4 

8.48e-5 

1.04e-5 

order 

- 

2.37 

3.01 

3.07 

3.03 

u  —  uh  HpO) 

1.48e-l 

3.03e-2 

4.02e-3 

5.27e-4 

6.88e-5 

order 

- 

2.29 

2.91 

2.93 

2.94 

1  u  —  u/i  L2(Q) 

5.23e-3 

6.40e-4 

5.59e-5 

4.63e-6 

3.27e-7 

order 

- 

3.03 

3.52 

3.59 

3.82 

|  Mr  (Ur)h  ||l2(0) 

0 

0 

0 

0 

0 

Me  —  (ue)h\  L2(Q) 

0 

0 

0 

0 

0 

\\p  -  PhWviSl) 

0 

0 

0 

0 

0 
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Table  15:  Annular  Poiseuillc  flow  convergence  rates  for  7  =  NURBS 


Polynomial  degree  k'  —  1 


h/h0 

1/2 

1/4 

1/8 

1/16 

1/32 

||u-Uh||V(h) 

1.73e-l 

9.88e-2 

4.91e-2 

2.35e-2 

1.13e-2 

order 

- 

0.81 

1.01 

1.06 

1.06 

u  —  Uh  hRh) 

8.20e-l 

5.29e-l 

2.90e-l 

1.49e-l 

7.46e-2 

order 

- 

0.65 

0.87 

0.96 

1.00 

u  —  u^l  L2(Q) 

4.89e-2 

2.12-2 

6.70e-3 

1.95e-3 

5.37e-4 

order 

- 

1.21 

1.66 

1.78 

1.86 

\\ur  —  (ur)h\  L2(f2) 

0 

0 

0 

0 

0 

\\ue  ~  (ue)h\ U2(fi) 

0 

0 

0 

0 

0 

\\p-Ph\\LHQ) 

0 

0 

0 

0 

0 

Polynomial  degree  k'  —  2 


h/ho 

1/2 

1/4 

1/8 

1/16 

1/32 

|  U  —  Uh||v(h) 

8.77e-2 

2.59e-2 

6.38e-3 

1.46e-3 

3.37e-4 

order 

- 

1.76 

2.02 

2.13 

2.12 

u  —  uh  HRO) 

4.34e-l 

1.39e-l 

3.63e-2 

8.86e-3 

2.16e-3 

order 

- 

1.64 

1.94 

2.03 

2.04 

u  —  uh  L2(Q) 

2.32e-2 

3.98e-3 

5.33e-4 

7. 14e-5 

9.46e-6 

order 

- 

2.54 

2.90 

2.90 

2.92 

||  "Uj.  (Ur)h\\L2(fl) 

0 

0 

0 

0 

0 

\We  ~  (ue)h\ \l2(Q) 

0 

0 

0 

0 

0 

\\p-Ph\\l?(Sl) 

0 

0 

0 

0 

0 

Polynomial  degree  k'  —  3 


h/ho 

1/2 

1/4 

1/8 

1/16 

1/32 

u  —  Uh||v(h) 

2.95e-2 

5.72e-3 

7.10e-4 

8.48e-5 

1.04e-5 

order 

- 

2.37 

3.01 

3.07 

3.03 

u  —  hRh) 

1.48e-l 

3.03e-2 

4.02e-3 

5.27e-4 

6.88e-5 

order 

- 

2.29 

2.91 

2.93 

2.94 

1  u  —  u/i  L2(Q) 

5.23e-3 

6.40e-4 

5.59e-5 

4.63e-6 

3.27e-7 

order 

- 

3.03 

3.52 

3.59 

3.82 

\\ur  iP'r)h  ||l2(0) 

0 

0 

0 

0 

0 

\\uq  —  {ue)h\  L2(Q) 

0 

0 

0 

0 

0 

\\P  ~  PhWmil) 

0 

0 

0 

0 

0 
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Figure  20:  Lid-driven  Stokes  flow  in  a  two-dimensional  cavity:  Problem  setup. 

9  Benchmark  Problems 

In  this  section,  we  investigate  the  effectiveness  of  our  methodology  as  applied  to  two 
standard  Stokes  benchmark  problems:  two-dimensional  and  three-dimensional  lid- 
driven  cavity  flow.  In  addition,  we  investigate  the  effectiveness  of  our  methodology  for 
a  Darcy-dominated  generalized  Stokes  flow  problem  characterized  by  sharp  boundary 
layers.  As  in  the  last  subsection,  we  choose  Nitsche’s  penalty  constant  as  Cpen  = 
5 (k'  +  1)  in  all  of  the  following  numerical  tests,  and  we  employ  uniform  parametric 
meshes,  linear  parametric  mappings,  and  B-spline  spaces  of  maximal  continuity. 

9.1  Two-Dimensional  Lid-Driven  Cavity  Flow 

Two-dimensional  lid-driven  cavity  flow  is  one  of  the  classical  verification  tests  for 
numerical  discretizations  of  incompressible  flow.  The  setup  for  this  flow  problem  is 
elaborated  in  Figure  20.  The  reaction  coefficient  a  and  the  applied  forcing  f  are  set 
to  be  zero.  The  left,  right,  and  bottom  sides  of  the  cavity  are  fixed  no-slip  walls  while 
the  top  side  of  the  cavity  is  a  wall  which  slides  to  the  right  with  velocity  magnitude 
U.  For  the  computations  here,  H  and  U  are  set  to  be  1.  The  pressure  and  stress  fields 
associated  with  this  flow  experience  corner  singularities  which  impede  the  convergence 
of  numerical  methods  and  expose  unstable  velocity/pressure  pairs.  In  fact,  the  exact 
velocity  solution  does  not  even  lie  in  H1(fi).  Instead,  it  lies  in  the  Sobolev  space 
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Wll9(fi)  where  1  <  q  <  2.  The  velocity  held  is  additionally  characterized  by  a 
primary  vortex  near  the  center  of  the  cavity  and  an  infinite  sequence  of  so-called 
Moffatt  eddies  of  decreasing  size  and  intensity  in  the  lower  left  and  right  corners  of  the 
cavity  [47] .  Of  primary  interest  is  how  well  our  discretization  procedure  approximates 
the  smooth  portions  of  the  how.  Of  secondary  interest  is  how  well  our  discretization 
procedure  resolves  the  corner  singularities.  We  do  not,  however,  employ  any  special 
techniques  such  as  singular  finite  elements  to  handle  the  singularities. 

We  have  computed  approximations  of  two-dimensional  lid-driven  cavity  how  using 
divergence-conforming  B-splinc  discretizations  of  varying  mesh  size  and  polynomial 
degrees  k'  =  1,2,3.  The  computed  streamlines  for  two  of  these  approximations  are 
presented  in  Figure  21.  The  streamlines  corresponding  to  the  approximation  defined 
on  the  fine  mesh  {k'  —  1,  h  —  1/128)  are  virtually  indistinguishable  from  well-accepted 
benchmark  solutions  [51].  The  streamlines  corresponding  to  the  approximation  de¬ 
fined  on  the  coarse  mesh  ( k '  —  1,  h  —  1/16)  closely  resemble  the  fine  mesh  streamlines 
in  the  interior  of  the  domain.  In  the  four  corners  of  the  domain,  the  approximation 
on  the  coarse  mesh  exhibits  visible  numerical  error  due  to  lack  of  resolution.  These 
results  indicate  that  our  methodology  suffers  from  minimal  pollution  error  [4].  It  is 
hypothesized  that  this  is  due  to  the  local  stability  and  approximation  properties  of 
B-splines.  To  further  highlight  how  well  the  approximation  on  the  coarse  mesh  ap¬ 
proximates  the  solution  in  the  interior  of  the  domain,  we  have  plotted  the  value  of  the 
first  component  of  the  velocity  held  along  the  vertical  center  line  in  Figure  22(a)  and 
the  value  of  the  second  component  along  the  horizontal  center  line  in  Figure  22(b) 
for  the  both  the  coarse  mesh  and  fine  mesh  solutions.  It  should  be  mentioned  that 
we  have  captured  the  first  Moffatt  eddy  in  the  two  lower  corners  with  both  meshes, 
and  we  have  observed  a  second  Moffatt  eddy  for  meshes  of  size  h  <  1/256. 

We  finish  our  discussion  of  the  two-dimensional  lid-driven  cavity  problem  by  an¬ 
alyzing  how  well  our  methodology  resolves  the  corner  singularities.  To  do  so,  we 
compare  the  results  of  our  methodology  with  the  highly  accurate  pseudospectral  re¬ 
sults  given  in  [10].  These  pseudospectral  results  were  obtained  using  a  subtraction  of 
the  leading  terms  of  the  asymptotic  solution  of  the  Stokes  equations  in  the  vicinity 
of  the  corners  in  order  to  exactly  represent  the  corner  singularities.  In  Table  16, 
we  compare  the  vorticity  (cu  =  curlu)  given  by  our  numerical  methodology  with  the 
pseudospectral  vorticity  near  the  upper  right  corner  of  the  cavity.  We  see  that  the 
vorticity  associated  with  our  numerical  methodology  slowly  converges  to  the  con¬ 
verged  pseudospectral  solution.  We  further  compare  our  computed  vorticities  with 
the  vorticity  obtained  with  a  highly- refined  finite  difference  solution  [32],  We  find 
that,  for  k'  —  1,  our  solutions  are  more  accurate  than  the  finite  difference  solution  for 
h  <  1/64,  and  for  k!  =  2,  our  solutions  are  more  accurate  than  the  finite  difference 
solution  for  h  <  1/32.  For  k'  =  3,  our  solutions  are  more  accurate  than  the  finite 
difference  solution  at  all  resolutions. 
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(a) 


(b) 


Figure  22:  Lid-driven  Stokes  flow  in  a  two-dimensional  cavity:  (a)  Value  of  the  first 
component  of  the  velocity  field  along  the  vertical  center  line  for  k'  —  1,  (b)  Value  of 
the  second  component  of  the  velocity  field  along  the  horizontal  center  line  for  k'  =  1. 
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Table  16:  Lid-driven  Stokes  flow  in  a  two-dimensional  cavity:  Convergence  of  vorticity 
at  the  point  (x  =  (1,0.95)). 


Polynomial  degree  k'  =  1 


Method 

UJ 

B-spline,  h  =  1/16 

-0.80995 

B-spline,  h  =  1/32 

14.34482 

B-spline,  h  =  1/64 

19.04468 

B-splinc,  h  =  1/128 

23.29179 

B-splinc,  h  =  1/256 

25.32238 

Pseudospectral  (Ref.  [10]) 

27.27901 

Finite  Difference  (Ref.  [32]) 

18.08 

Polynomial  degree  k'  - 

=  2 

Method 

UJ 

B-spline,  h  =  1/16 

11.06384 

B-spline,  h  =  1/32 

31.81761 

B-spline,  h  =  1/64 

32.81972 

B-splinc,  h  =  1/128 

26.48645 

B-splinc,  h  =  1/256 

27.34395 

Pseudospectral  (Ref.  [10]) 

27.27901 

Finite  Difference  (Ref.  [32]) 

18.08 

Polynomial  degree  k'  - 

=  3 

Method 

UJ 

B-spline,  h  =  1/16 

29.86220 

B-spline,  h  =  1/32 

25.00897 

B-spline,  h  =  1/64 

29.92944 

B-splinc,  h  =  1/128 

28.75895 

B-splinc,  h  =  1/256 

27.52637 

Pseudospectral  (Ref.  [10]) 

27.27901 

Finite  Difference  (Ref.  [32]) 

18.08 
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Figure  23:  Lid-driven  Stokes  flow  in  a  three-dimensional  cavity:  Problem  setup. 

9.2  Three-Dimensional  Lid-Driven  Cavity  Flow 

While  three-dimensional  lid-driven  cavity  flow  is  encountered  much  less  often  in  the 
literature  than  its  two-dimensional  counterpart,  we  believe  it  is  still  an  interesting  test 
case  for  numerical  discretizations  of  incompressible  flow.  In  the  three-dimensional  set¬ 
ting,  cavity  flow  is  characterized  by  the  presence  of  both  edge  and  corner  singularities. 
The  problem  setup  is  illustrated  in  Figure  23.  Again,  for  the  computations  here,  o 
and  f  are  set  to  be  zero.  Every  side  of  the  cavity  except  the  top  is  assumed  to  be 
a  stationary  no-slip  wall,  and  the  top  side  of  the  cavity  is  assumed  to  be  a  wall 
which  slips  to  the  right  with  velocity  magnitude  U .  The  streamlines  resulting  from  a 
simulation  of  three-dimensional  lid-driven  cavity  flow  using  divergence-conforming  B- 
splines  of  degree  k!  =  2  on  a  mesh  with  32  x  32  elements  are  illustrated  in  Figure  24. 
Note  that  the  three-dimensional  streamlines  resemble  the  two-dimensional  stream¬ 
lines  along  the  slice  y/H  =  1/2.  To  examine  how  well  our  discretization  technique 
performs  on  coarse  meshes,  we  have  compared  the  centerline  values  of  our  velocity 
field  along  the  slice  y/H  =  1/2  for  both  a  coarse  and  fine  mesh  approximation  in 
Figure  25.  From  the  figure,  we  see  that  the  coarse  and  fine  mesh  centerline  velocity 
fields  are  nearly  indistinguishable.  This  indicates  that  our  methodology  suffers  from 
minimal  pollution  error. 
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Figure  24:  Lid-driven  Stokes  flow  in  a  three-dimensional  cavity:  Computed  flow 
velocity  streamlines  for  k'  —  2  and  h  —  1/32  colored  by  velocity  magnitude. 


9.3  Darcy-Dominated  Flow  with  Sharp  Boundary  Layers 

As  a  final  numerical  example,  we  consider  a  reaction-dominated  generalized  Stokes 
problem  subject  to  sharp  boundary  layers.  The  problem  is  posed  on  the  unit  square. 
Homogeneous  no-slip  and  no-penetration  boundary  conditions  are  enforced  along  the 
boundary  of  the  square,  and  an  external  forcing  of  the  form 


f  = 


sin(7nr)  cos(7 ry) 
—  cos(7nr)  sin(7n/) 


is  applied.  The  reaction  term  a  is  set  to  be  equal  to  1  while  the  viscosity  v  is  set  to 
be  equal  to  10“6.  With  these  choices,  the  resulting  flow  field  has  a  vortical  structure 
in  the  interior  of  the  domain  and  sharp  boundary  layers  along  the  entire  boundary  of 
the  domain. 

We  have  simulated  this  flow  problem  using  both  our  discretization  technique  as 
well  as  the  Q2/Qi  Taylor-Hood  velocity/pressure  pair  [36].  The  Taylor- Hood  veloc¬ 
ity/pressure  pair  is  one  of  the  most  popular  finite  elements  for  generalized  Stokes  flow, 
but  it  is  known  not  to  be  robust  in  the  Darcy  limit.  Specifically,  the  accuracy  of  the 
Q2/Q1  Taylor-Hood  velocity  field  reduces  to  first-order  in  the  Darcy  limit  [33].  This  is 
due  to  a  lack  of  strong  coercivity  (in  the  kernel)  with  respect  to  the  H(div;  Q)  -norm. 
We  have  visualized  the  computed  velocity  vectors  corresponding  to  divergence-free 
B-splines  of  degree  k'  =  1  and  the  Q2/Qi  Taylor-Hood  velocity-pressure  pair  on  a 
16  x  16  element  mesh  in  Figures  26(a)  and  (b).  These  two  quiver  plots  seem  to  indi¬ 
cate  that  the  two  discrete  flow  fields  are  nearly  identical  except  in  a  small  region  near 
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x/H 

(b) 

Figure  25:  Lid-driven  Stokes  flow  in  a  three-dimensional  cavity:  (a)  Value  of  the  first 
component  of  the  velocity  field  along  the  line  (x/H,y/H)  =  (1/2, 1/2)  for  k'  =  1, 
(b)  Value  of  the  third  component  of  the  velocity  field  along  the  line  (y / H,  z/ H)  = 
(1/2, 1/2)  for  k'  =  1. 
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the  boundary  of  the  domain.  However,  when  we  view  only  the  first-component  of  the 
velocity  held  as  in  Figures  27  and  28,  we  see  that  the  divergence- free  B-spline  solu¬ 
tion  is  monotone  while  the  Taylor-Hood  solution  suffers  from  spurious  oscillations. 
Finally,  we  have  visualized  the  divergence  of  the  discrete  how  held  corresponding  to 
the  Taylor-Hood  solution  in  Figure  29.  From  this  figure,  we  see  that  the  Taylor-Hood 
solution  is  characterized  by  strong  expansion  and  compression  in  the  four  corners  of 
the  domain.  When  coupled  with  a  transport  solver,  this  ultimately  leads  to  methods 
with  unphysical  species  production  [46].  On  the  other  hand,  the  discrete  how  held 
corresponding  to  the  B-spline  solution  is  pointwise  divergence-free. 

10  Conclusions 

In  this  paper,  new  divergence-conforming  B-spline  discretizations  have  been  presented 
for  the  generalized  Stokes  equations.  A  collection  of  stability  and  error  estimates  have 
been  derived  which  are  robust  with  respect  to  the  parameters  of  generalized  Stokes 
how,  and  these  theoretical  results  have  been  confirmed  and  supplemented  by  numer¬ 
ical  simulations  of  problems  with  known  analytical  solutions.  These  discretizations 
have  been  also  applied  to  the  simulation  of  a  number  of  benchmark  problems  where 
the  advantages  of  the  new  methodology  over  classical  methods  have  been  highlighted. 
In  future  work,  we  plan  to  extend  the  present  developments  to  the  Navier-Stokes 
equations  and  generalize  to  locally-refined  meshes. 
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Figure  26:  Darcy-dominated  flow  with  boundary  layers:  Computed  velocity  fields:  (a) 
Divergence-free  B-splines  of  degree  k!  =  1,  (b)  Q2/Qi  Taylor-Hood  velocity/pressure 
pair. 
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Figure  27:  Darcy-dominated  flow  with  boundary  layers:  First-component  of  the  ve¬ 
locity  field:  (a)  Divergence-free  B-splines  of  degree  k!  =  1,  (b)  Q2/Qi  Taylor-Hood 
velocity/pressure  pair. 
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Figure  28:  Darcy-dominated  flow  with  boundary  layers:  Cross-section  view  of  the 
first  velocity  component  from  the  y  =  1  line:  (a)  Divergence-free  B-splines  of  degree 
k'  =  1,  (b)  Q2/Qi  Taylor-Hood  velocity/pressure  pair. 
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Figure  29:  Darcy-dominated  flow  with  boundary  layers:  Divergence  of  the  discrete 
flow  field  corresponding  to  the  Q2/Qi  Taylor-Hood  velocity /pressure  pair. 
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